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ABSTRACT 

Two  GLR  algorithms  are  developed  for  detecting  and  estimating 
abrupt  maneuvers  (i.e.,  jumps  in  control  values)  in  discrete 
linear  stochastic  systems.  A  jump  error  state  variable  concept 
is  used  to  derive  a  decomposition  of  the  conventional  maneuver 
signature  matrix  into  a  new  maneuver  signature  matrix  which 
is  independent  of  the  jump  time  and  another  matrix  which  depends 
only  on  the  jump  time.  The  product  of  the  latter  matrix  with 
the  jump  vector  is  shown  to  provide  a  constant  jump  error  state 
variable.  The  nondependency  of  the  new  maneuver  signature 
matrix  on  the  jump  time  and  the  transformation  of  the  jump 
error  state  to  a  cone  provides  for  the  development  of  GLR 
algorithms  with  considerably  reduced  computational  and  storage 
requirements.  The  new  algorithms  avoid  much  of  the  updating 
and  storage  of  large  matrices  for  past  observation  times. 

President,  Practical  Sciences, Inc. ,  40  Long  Ridge  Road, 

Carlisle,  MA  01741 


In  particular,  the  multiplications  requirement  is  reduced  to 
3  2 

the  order  of  n  +  3Mn  m  for  general  linear  stochastic  systems 
2 

and  Mn  m  for  such  systems  without  process  noise.  Here,  n  is 
the  dimension  of  the  state,  m  the  dimension  of  the  measurement 
vector  and  M  is  number  of  candidate  jump  times  in  the  past. 

The  two  algorithms  have  practical  application  to  the  engage¬ 
ment  problem  between  an  anti-ship  cruise  missile  and  a  ship 
defense  interceptor.  The  output  of  the  algorithms  provides 
information  on  which  the  players  in  a  differential  game  of 
partial  information  and  noisy  observation  may  base  and  design 
optimal  strategies  for  maximizing  payoff  functions  such  as 
survivability  and  kill. 
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1.  BACKGROUND  AND  INTRODUCTION 


The  engagement  between  an  anti-ship  cruise  missile  (ASM) 
and  a  ship  defense  interceptor  is  a  differential  game  of 
partial  information  and  noisy  observations.  The  ASM  represents 
the  offense  (the  evader)  and  the  interceptor  represents  the 
defense  (the  pursuer).  The  ASM  desiring  to  optimize  its 
survivability  while  enroute  to  its  ship  target  employs  endgame 
maneuvers  for  the  purpose  of  evading  the  ship' s  defense  interceptor 
guided  by  active  and/or  passive  radar.  The  ASMs  are  highly 
sophisticated  and  maneuvering  missiles  that  have  speeds  ranging 
from  subsonic  to  supersonic,  altitudes  ranging  from  sea  skimmer 
to  upper  limits  of  winged  aircraft,  and  guidance  systems 
ranging  from  simple  line-of -sight  to  more  complex  combinations 
of  passive  and  active  electronic  systems.  Their  maneuvering 
characteristics  are  governed  by  airframe,  propulsion  and 
guidance  parameters  —  drag,  lift,  thrust  and  four  guidance 
parameters,  [1]  and  [2].  The  thrust  and  guidance  parameters 
are  like  step  fucntions  in  time  (i.e.,  piecewise  constant 
controls) .  They  take  on  one  value  during  one  portion  of 
the  ASM's  trjectory  and  jump  to  other  values  during  other 
portions  of  the  trajectory.  For  example,  the  guidance  parameters 
jump  in  value  at  the  start  of  pop-ups,  dives, pull-outs,  turn-downs 
and  turn-on  of  seeker  homing.  The  thrust  parameter  of  some 
ASMs  jumps  (i.e.,  drops)  in  value  during  high  altitude  dives. 
Consequently,  the  dynamics  of  an  ASM  are  representable  as 
a  stochastic  system  governed  by  piecewise  constant  controls. 

In  linear  discrete  form  the  dynamics  are  modeled  as: 


(1) 


System  Dynamics  of  ASM* 

X  (k+1)  =  A(k+l,k)  X(k)  +  B (k+l,k)  (u(k)  +  Auq  6qk)  +  r(k)w(k) 
u(k+l)  *  u(k)  +  6qk 


where  X  is  the  state  vector,  u  is  the  control,  A  uq  is  the 
jump  in  control  at  time  q,  6qk  is  the  Kronecker  delta,  and  T 
is  the  system  noise  coefficient  matrix.  The  matrix  A  is 
the  state  transition  matrix  and  B  is  the  input 
matrix  for  the  control.  The  jump  time  q  and  the  jump  magnitude 
Auq  are  unknown  to  the  defense  interceptor.  Using  active 
and/or  passive  radar  sensors  the  defense  interceptor  measures 
a  partial  state  of  the  ASM.  These  observations  are  modeled 
as: 


Sensor  Equation  of  Interceptor 

z  (k)  =  h  (k)  X  (k)  +  v(k)  (3) 

where  z  is  the  measurement  vector  and  h  is  the  measurement 
matrix.  The  noise  sequences  w  and  v  are  zero-mean,  independent, 
white  Gaussian  sequences  with  covariances  defined  by 


♦Actually,  there  are  multiple  jumps  Auq.at  Limes  qj.;  this 
is  easily  indicated  by  the  substitution1 of  \  Au<Ii  6<U* 

into  (1)  and  (2).  Equations  (1)  and  (2)  represent  each 

jump  in  turn. 


2 


E  (w(k)  wT  (j)}  *  Q  (k)  Sfcj 


E  {0  00  vT(j)}  e  R  <k)  6  k  j 

where  e{#>  denotes  the  expectation  and  the  matrix  R(k)  is 
bounded  positive  definite.  The  inital  state  X(o)  is  normally 

a 

distibuted  with  mean  X(o)  and  covariance  P(0);  we  make  a 
similar  assumption  for  u(0). 

The  system  dynamics  of  the  interceptor  is  of  a  form 
similar  to  (1)  and  (2).  In  general,  the  sensor  equation 
of  an  ASM  may  be  considered  to  be  of  the  form  (3)  but  currently, 
in  practice,  the  ASM  has  no  sensor  with  which  to  observe 
the  location  or  presence  of  an  interceptor  (i.e.,  it  is  blind 
to  approaching  interceptors)  even  though  it  has  the  sensors 
for  acquiring  and  homing  in  on  its  ship  target.  In  such 
a  case  the  measurement  matrix  of  the  ASM  is  zero.  There 
are  other  engagement  scenarios  in  which  the  offense  (e.g., 
bomber)  is  not  blind  to  approaching  interceptors.  Therefore, 
in  general,  we  are  interested  in  the  class  of  differential 
games  of  partial  information  and  noisy  observations  in  which 
the  dynamics  of  a  player's  system  are  modeled  by  (1)  and 
(2)  and  the  observation  equation  of  its  opponent  is  modeled 
by  (3). 
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In  the  above  engagement  the  ASM  desires  to  maximize 
its  survivability  using  evasive  maneuvers  while  enroute  to 
its  ship  target.  The  interceptor  desires  to  maximize  its 
probability  of  kill.  In  its  pursuit  of  optimality  the  interceptor 
is  faced  with  two  tasks.  The  first  is  the  development  of 

/V  A 

an  estimator  for  obtaining  the  optimal  estimates  X  and  u 
of  the  state  X  and  the  control  u.  Secondly,  it  has  the  task 
of  deriving  optimal  strategies  based  on  the  optimal  estimates. 

In  this  paper  we  address  the  first  task  and  present  a  filtering 
algorithm  for  obtaining  the  optimal  estimates  of  X  and  u. 

This  estimation  problem  is  that  of  detecting  and  estimating 
abrupt  changes  in  stochastic  dynamical  systems.  A  survey 
of  estimation  methods  is  given  in  Willsky  [3]. 

One  of  the  most  attractive  and  promising  methods  for 
detecting  and  estimating  jumps  in  linear  stochastic  systems 
is  the  generalized  likelihood  ratio  (GLR)  method,  [4] -[6], 

The  GLR  method  processes  the  residuals  from  a  Kalman  filter 
and  computes  the  maximum  likelihood  estimates  of  the  jump 
time  and  the  jump  magnitude.  Using  these  estimates  it  evaluates 
the  log-likelihood  ratio  for  jump  versus  no  jump.  A  jump 
is  declared  if  the  evaluated  ratio  is  larger  than  a  set 
threshold.  The  implementation  of  the  GLR  requires  a  linearly 
growing  bank  of  matched  filters  in  order  to  compute  the  maximum 
likelihood  estimate  of  the  jump  time,  [3}  and  [7] .  A  recursive 
GLR  algorithm,  [8]  and  [9] ,  has  been  developed  that  reduces 
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the  computational  burden  by  reducing  or  avoiding  the  requirement 
for  matrix  inversion  (in  computing  the  jump  magnitude  and 
evaluating  the  log-likelihood  ratio)  at  each  possible  jump 
time  in  the  past.  This  reduction  was  obtained  by  modifying 
the  GLR  algorithm  so  that  the  covariance  of  the  predicted 
measurement  residual  is  to  be  inverted  rather  than  the  information 
matrix  of  the  jump  variable.  The  largest  dimension  of  the 
matrix  to  be  inverted  is  at  most  equal  to  the  dimension  of 
the  correlated  components  of  the  predicted  measurement  residual, 

[10]  and  [11]. 

In  its  latest  stage  of  development  the  GLR  method  still 
requires  at  each  new  observation  time  the  computation  and 
storage  of  several  matrices  for  each  observation  time  in 
the  past.  Herein,  we  derive  two  GLR  algorithms  I  and  II  which 
are  based  on  a  decomposition  of  the  failure  signature 
matrix  of  [5]  which  we  term  herein  the  maneuver  signature  matrix. 

The  decomposition  provides  a  new  maneuver  signature  matrix 
which  is  independent  of  the  jump  time  (maneuver  start  time) . 

This  nondependency  reduces  considerably  the  storage  and  computational 
requirements  of  the  GLR  method  by  reducing  the  requirement 
to  update  and  restore  large  matrices  at  each  current  observation 
time  for  past  observation  times.  The  computational  burden 
of  previous  GLR  algorithms  is  to  a  great  extent  the  direct 
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dependency  of  the  maneuver  signature  matrix  on  the  jump  time. 

Our  algorithms  I  and  II  avoid  much  of  the  updating  and  storage 
of  large  matrices  for  past  observation  times. 

The  GLR  algorithms  I  and  II  are  developed  using  the  concept 
of  a  jump  error  A-state  variable  Ax.  It  is  introduced  as  the 
difference  between  a  "jump  is  known"  filter  estimate  and  a 
"jump  free"  filter  estimate  of  the  state.  The  evolution  of 
Ax,  after  a  jump,  is  governed  by  a  linear  state  equation.  The 
residuals  of  the  "jump  free"  filter  provides  the  noisy  linear 
measurements  of  Ax.  The  variable  Ax  is  non-singular ly  transformed 
into  a  new  A-state  variable  Ay  which  is  constant.  It  is  that 
transformation  which  provides  the  new  maneuver  signature  matrix 
which  is  independent  of  the  jump  time. 


Our  work  is  built  on  that  of  Friedland's  bias  filtering 
technique  [20],  Willsky  and  Jones'  GLR  technique  [6]  and  Chang 
and  Dunn's  recursive  GLR  algorithm  [9], 
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2.  DISCRETE  LINEAR  STOCHASTIC  MODEL  WITH  STATE  JUMP: 

A  POSTERIORI  JUMP 

For  ease  of  presentation  and  generality  we  augment  the 
state  X(k)  with  the  control  u(k)  of  (1)  and  (2)  and  define 
the  variable  x(k)  as  the  augmented  state  vector  composed 
of  X(k)  and  u(k).  In  this  case  Equations  (1),  (2),  and  (3) 
become 


x(k+l)  =  9  (k+l,k)  (x  (k)  +  AXq  6qk)  +  T  (k)w(k)  (4) 

z  (k)  =  H  (k)  x  (k)  +  v  (k)  (5) 

where 

A(k+l,k)  B (k+l,k) 

$(k+l,k)  =  [  ]  (6) 

0  I 

A  xq  =  D  A  uq  (7) 

0 

D  =  [  ]  (8) 

I 

H (k)  =  [h (k)  0]  (9) 


and  where  T  (k)  is  redefined  as  the  augmented  matrix 
T  (k) 

(  1  (10) 
0 

A  description  of  the  variables  and  their  dimensions  are  given 
in  Table  1.  We  assume  that  the  linear  system  (4)  and  (5) 
is  observable. 
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The  a  posteriori  jump  is  used  in  the  formulation  of 
(1),  (2),  and  (3).  The  jump  A  Xq  occurs  at  time  q  but  it 
does  not  appear  in  the  measurement  until  time  q+1;  that  is, 
the  jump  occurs  right  after  the  measurement  z(q). 

The  optimal  state  estimator  for  a  discrete  linear  stochastic 
system  without  jump  (Axq  *  0)  is  given  by  the  discrete  Kalman-Bucy 
filter,  [ 12 ] — [ 15 ] -  The  equations  of  the  filter,  [16],  are 
given  in  Table  2  for  reference  purposes.  The  filter  variables 
are  described  in  Table  3. 

In  this  report  we  treat  the  general  problem  as  defined  by 
Eqs.  (4)  and  (5)  in  which  the  unknowns  to  be  estimated  and  detected 
are  Ax  and  q.  That  is,  we  take  the  quantities  $  ,  Ax  ,  H 

,  .  i  4  'a 

and  r  to  be  arbitrary,  not  necessarily  satisfying  (6),  (7),  (9) 
and  (10)  . 
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TABLE  1 


SYSTEM  VARIABLES  FOR  THE  AUGMENTED  SYSTEM 


VARIABLE _ DEFINITION _ DIMENSION 


x(k) 

State  vector 

nxl 

$(k+l,k) 

State  transition  matrix  from 
k  to  k+1 

nxn 

T(k) 

Q(k) 

A  Xq 

System  noise  coefficient  matrix 
System  noise  covariance  matrix 
Jump  in  state  at  time  q 

nxr 

rxr 

nxl 

z(k) 

Measurement  at  time  k 

mxl 

H  (k) 

Measurement  matrix 

mxn 

R(k) 

Measurement  noise  covariance 
matrix 

mxm 

w(k) 

Gaussian  white  system  noise 

rxl 

v(k) 

Gaussian  white  measurement  noise 

mxl 

D 

Jump  coefficient  matrix 

nxp 

A  uq 

Jump  in  control  at  time  q 

pxl 

u(k) 

Control  vector 

pxl 
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TABLE  2 


DISCRETE  KALMAN-BUCY  FILTER  EQUATIONS* 


x (k+1 1 k)  =  $(k+l,k)  x{k)  (i) 

P  (k+1  Ik)  =  *(k+l,k)  P  (k)  $T(k+i,k)  +r(k)  Q(k)rT(k)  (ii) 

y  (k)  =  z  (k)  -  H(k)  x(klk-l)  (iii) 

V(k)  *  H (k) P (k | k-1)  HT(k)  +  R(k)  (iv) 

K(k)  *  P (k 1 k-1)  HT(k)  V-itk)  (v) 

/s  /v 

x  (k)  *  x.(k|k-l)  +  K(k)  Y(k)  (vi) 

P ( k )  ■  (I  -  K (k)  H (k ) ]  P (k | k-1)  (vii) 


*The  usual  notations  x  (k I k)  and  P(k|k)  have  been  shortened 

A  A 

to  x  (k)  and  P(k).  The  random  variable  x(k|j)  is  the  optimal 
estimate  of  x (k)  based  on  all  the  measurements  Z(j)  ={z(l)f 
z(2),...,  z(j)}  .  The  superscript  "T"  denotes  transpose 
and  "-1"  denotes  inverse.  The  identity  matrix  is  denoted 
by  I. 
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TABLE  3 


FILTER  VARIABLES 


VARIABLE  DEFINITION _ DIMENSION 


A 


x(k) 

State  estimate  at  k  given  Z(k) 

nxl 

P(k) 

Covariance  matrix  of  the  error 
in  x  (k) 

nxn 

*  ik+ll 

Ik) 

State  estimate  at  k+1  given  Z(k) 

nxl 

P  (k+1 1 

Ik) 

Covariance  matrix  of  the  error 
in  x  (k+1 Ik) 

nxn 

Y(k) 

Predicted  measurement  residual 

mxl 

V(k) 

Covariance  of  y  (k) 

raxm 

K(k) 

Filter  (Kalman)  gain  matrix  at  k 

nxm 

II 


3.  FORMULATION  OF  THE  A-SYSTEM:  JUMP  ERROR  STATE  EQUATION 


Consider  the  following  filtering  conditions  for  the  Kalman- 
Bucy  filter: 

:  There  is  no  jump  in  state  and  no  jump  in  assumed  by 
the  filter. 

H^:  There  is  a  jump  in  state  but  the  filter  is  unaware 
that  a  jump  has  taken  place  and  it  operates  as  if 
the  jump  is  zero.  The  filter  is  referred  to  as  the 
"jump  free"  filter. 

H2:  There  is  a  jump  in  state,  the  jump  is  known  to  the 
filter  and  the  jump  information  (time  and  magnitude) 
is  made  use  of  in  the  filter.  The  filter  is  called 
the  "jump  is  known"  filter. 


The  Kalman-Bucy  filter  is  optimal  for  Conditions  ^ 
and  H2  and  nonoptimal  for  H^.  Condition  H2  is  ideal  but 
does  not  occur  in  practice.  Condition  is  the  real  world 
condition.  We  are  faced  with  the  problem  of  accounting  for 
the  jump  after  it  has  occurred  and  has  been  detected.  On 
the  one  hand  we  do  not  wish  to  degrade  the  optimal  performance 
of  the  Kalman-Bucy  filter  by  operating  it  in  an  "after-jump" 
mode  when  no  jump  has  occurred.  On  the  other  hand  we  desire 
optimal  estimates  of  the  state  after  the  jump  has  occurred. 

Because  of  the  delay  between  the  occurrence  of  the  jump  and 
its  detection  the  output  of  the  filter  is  nonoptimal  during 
this  delay.  Since  we  have  in  practice  this  sequence  of  nonoptimal 
estimates  we  would  like  to  compensate  it  with  an  additional 
estimate  and  obtain  an  optimal  estimate.  This  is  precisely 
the  property  of  the  a -system.  The  optimal  estimate  of  the 
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A-state  added  to  the  Kalman-Bucy  filter  estimate  under  condition 
is  an  optimal  estimate  of  the  state  beyond  the  jump. 

Let  x i (k)  represent  the  state  for  the  case  that  the 
jump  magnitude  Axq  is  zero  (i.e.,  there  is  no  jump).  Let 
x_2  (k)  represent  the  state  for  the  case  that  there  is  a  jump 

A  M  A* 

AXg  and  its  magnitude  may  be  nonzero.  Let  Xj»  and  *2. 

denote  the  Kalman-Bucy  filter  estimates  of  the  state  for 
the  Conditions  and  H2,  respectively.  The  jelationships 

between  these  estimates  are  depicted  in  Figure  1  for  an  anti¬ 
ship  missile  defense  scenario  in  which  the  ASM  employs  an 
endgame  pop-up  maneuver.  Under  Condition  the  Kalman- 
Bucy  filter  estimate  optimally  tracks  •  Under  Condition 
H2  the  Kalman-Bucy  filter  estimate  x-2  optimally  tracks  x:2. 

But  under  condition  the  nonoptimal  estimate  tracks 
a  trajectory  between  Xj^  and  x2»  Since  the  covariance  of 
state  estimate  and  the  Kalman  gain  are  independent  of  the 
measurements  it  follows  that  the  gains  Kj  (k)  ,  K^  (k)  and  K2  (k) 
are  identical  and  that  the  covariances  Plf  and  P2  are 

A* 

identical  for  the  three  filtering  Conditions  H^,  and  H2. 

We  define  the  following  errors  of  (k) : 

A x(k)  *  *2  (k)  -  (k),all  k  (10] 

A w(k)  -  Xl(k)  -  ^(k),  all  k  (11] 

The  sum  of  the  two  errors  satisfy 
Ax  (k)  +Aw(k)  «  $(k,q)  fr  ,  k>q  U2] 
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ALTITUDE  (FT) 


KALMAN-BUCY  FILTER  ESTIMATES 


X^  -  NO  MANEUVER  *l  *  KALMAN  GAIN 

x2  -  MANEUVER  IS  KNOWN  H  =  MEASUREMENT  MATRIX 

X  -  MANEUVER  IS  UNKNOWN  *  =  STATE  TRANSITION  MATRIX 


4  3  2  1  0 

RANGE  TO  GQ.  (NMI ) 


FIGURE  1.  ASM  DEFENSE  SCENARIO:  JUMP  ERROR  STATE  Ax  SATISFIES 
LINEAR  EQUATION 

Ax (k+1)  =  [I-Ki (k+1)  H (k+1) ]  $(k+l,k)  Ax(k) 
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We  call  the  random  variable  ax  the  jump  error  state  or 
A-state  variable.  It  is  the  difference  between  the  Kalman- 
Bucy  filter  estimates  under  the  conditions  of  known  jump 
and  unknown  jump.  Since  the  addition  of  Ax(k)  and  x^fk)  give 

A-  A 

x2 (k)  we  would  expect  that  the  optimal  estimate  Ax(k)  added 

~  a, 

to  x^ (k)  gives  the  optimal  estimate  x(k)  of  (4)  and  (5): 


x(k)  «  x1  (k)  +  Ax(k) 


It  is  easy  to  show  that  ax  satisfies  a  linear  state 
equation  [8] .  The  initial  condition  for  &x  is  at  time  q: 

Ax(q)  *  a Xg 

^  ***' 

Axg  *  x2(q)  -  x^q) 

Let  k  >q.  From  Equations  (i)  ,  (iii) ,  and  (vi)  of  Table  2 

^  A 

we  observe  that  the  estimates  x^ (k)  and  x2 (k)  satisfy 
Xjlk)  *  A# (k,k-l)  x^k-1)  +  Kx  (k)  z(k) 


x2(k)  =a*  (k,k- 1)  x2(k-l)  +  K1(k)z(k) 
where 

A*  (k,k-l)  =  tI-Kx(k)  H (k)  3  *(k,k-l) 


Consequently,  subtracting  (16)  from  (17)  gives  the  linear 
equation 


Ax(k)  Ax(k-l) 


15 


with  initial  Condition  (14)  .  The  measurement  equation  for 
Ax  is  also  easily  derived,  under  Condition^,  the  a  posterior 
measurement  residual  az  00  is  given  by 

Az(k)  -  z(k)  -  H (k)  xx(k)  (20) 

A 

Adding  and  subtracting  the  term  H (k)  x2(k)  we  obtain  the 
measurement  equation  for  a* 

Az(k)  =  AH(k)  Ax(k)  +Av  00  (21) 

where 

A  H  (k )  =  H  (k)  (22) 

AvOO  -  z(k)  -  H (k)  x2(k)  (23) 

is  the,a  posteriori  measurement  residual  under  Condition  H2- 
Consequently,  AvOO  is  a  zero-mean  white  Gaussian  sequence 
with  covariance  defined  by 

E{Av(k)AvT(j)}  -  AR(k)6k;j  (24) 

where 

AR(k)  -  R(k)  V1“1(k)  R(k)  (25) 

recalling  that  and  V2  are  identical. 


Equations  (19)  and  (21)  constitute  the  linear  equations 
that  govern  the  error  Ax.  If  the  jump  time  q  wer'»  known 
and  the  initial  state  Ax(q)  were  normally  distributed  with 

/s 

mean  Ax(q)  and  covariance  aP(<2)  then  a  Kalman-Bucy  filter 
could  be  employed  to  estimate  Ax(k).  The  filtering  equations 
are  as  given  in  Table  4  for  such  a  case.  For  this  case  define 
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TABLE  4 


FILTERING  EQUATIONS  FOR  THE  A-SYSTEM 


Ax(k+l|k)  *  A*  (k+l,k)  Ax(k),  k>q  (i) 

AP(k+l|  k)  *  A$(k+l,k)  AP(k)  A*T(k+l,k),  k>  q  (ii) 

a 

Ay  (k)  =  Az(k)  -  AH(k)  Ax(klk-l),  k>q  (iii) 

AV(k)  *  AH (k)  AP  (k |  k-1)  AHT(k)  +  AR(k)  ,  k  >q  (iv) 

AK(k)  -  AP  (k|  k-1)  AHT  (k)  AV_1  (k)  ,  k>q  (v) 

A  A 

Ax(k)  =  Ax(k|k-1)+  AK(k)  Ay  (k)  ,  k>  q  (vi) 

AP(k)  -  [I  -  AK  (k)  AH  (k)  0  AP(k|k-l),  k>  q  (vii) 

where 

A«k,k-1)  *  [I  -  Kx(k)  H (k)  ]  *  (k,k-l)  (viii) 

AH  (k)  =  H  (k)  (ix) 

AR(k)  -  R(k)  V1’1(k)  R(k)  (x) 

Az(k)  *  z(k)  -  H  (k)  x^k)  (xi ) 
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(26) 


x  (q)  *  xi  (q)  +  Ax(q) 

P(q)  «  Pl(q)  +  AP(q)  (27) 


and  employ  a  Kalman-Bucy  filter  to  estimate  x(k),  k>q,  governed 

a 

by  (4)  and  (5).  We  know  that  the  resulting  estimates  x(k)  and 
P(k)  and  the  gain  K(k)  are  optimal.  On  the  other  hand,  we  can 
employ  the  filter  under  Condition  Hi  to  generate  the  estimates 
xi(k)  and  Pi(k),  k>q  and  we  can  employ  the  filter  of  Table  4 

A 

to  generate  the  estimates  Ax(k)  and  AP(k).  The  two  approaches 
are  equivalent.  That  is,  (this  is  shown  in  Appendices  C  and  D) , 


a  ~  A 


x  (k)  =  x1  (k)  +  Ax  (k) 

(28) 

P  (k)  =  P1(k)  +  AP(k) 

(29) 

For  such  a  case  the  optimal  estimate  of  x(k)  is  obtained 
by  summing  the  two  estimates  x^ (k)  and  Ax(k).  The  gains 
are  related  by  the  expression  (this  is  shown  in  Appendix  A) . 


K(k)  *  Kx(k)  +  AK  (k)  [I  -  H'(k)Ki  ( k>  ] 

(30) 

one  can  show  that 

E{A  e(k)  e^tk)  >  -  0,  k>  q,  i  =  1,2 

(31) 

where,  for  kj>q» 

e^ (k)  «  xA(k)  -  xi(k),  i*l,2 

(32) 

Ae(k)  *  Ax(k)  -  Ax(k) 

(33) 
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Note  that 


e1  (k)  =  e2(k),  x2 (k)  =  x{k)  and 
e(k)  =  e2(k)  +  Ae(k),  k  ^  q 
where 

e(k)  =  x(k)  -  x(k) 

From  (28)  we  note  that 
x  (k)  =  x2  (k)  +  Ax(k)  -Ax(k),  k  ^  q 


which  shows  that  the  optimal  estimate  x{k)  is  as  close  to 

A 

the  ideal  estimate  x2 (k)  (i.e.,  jump  is  known)  as  the  optimal 

/v 

estimate  Ax (k)  is  to  Ax(k). 

It  is  shown  in  Appendix  B  that  the  Predicted  Measurement 
Residual  covariances  are  related  by  the  expression 

V(k)  =  V1  (k)  R-1  (k)  AV(k)  R- 1  ( k )  Vx(k) 


(34) 

(35) 

(36) 


I 

I 

I  4.  TRANSFORMATION  OF  A  -STATE:  CONSTANT  £ -STATE 

! 

Define  the  nxn  matrix  \p(k)  ,  all  k,  as 

,|,(0)  =  I  (37) 

4,(k)  =  a$  (k,k-l)  4,(k-l),  k>  0  (38) 

The  matrix  ,|,(k)  is  positive  definite  for  all  k;  it  has  inverse 
4»-1(k).  The  &-state  equation  (19)  can  be  rewritten  as 

fix(k)  =  4»(k)  ^-1(k-l)  ^x(k-l),  k>  q  (39) 

we  define  a  new  A~state  variable  ay  as 

Ax(k)  =  ip  (k)  /\y(k),  k >  q  (40) 

This  transformation  of  the  ^-state  from  ^x  to  ay  results 
in  the  constant  state  equations 

Ay(k)  =  AY(k-l)  k>  q  (41) 

The  ^-measurement  matrix  ^H(k)  in  (22)  car.  l>e  rec-ifined  as 

AH(k)  =  H(k)  4,  (k)  (42) 

so  that  the  Equation (21)  in  terms  of  the  new  A-state  Ay  becomes 


Az  (k)  =  AH  (k)  Ay(k)  +Av  (k) 


(43) 


r'  ■Trjfclif 


Equations  (41)  and (43)  define  the  linear  A-system  for  the 
state  Ay.  Equation  (40)  gives  the  transformation  back  to 
the  ax  state. 

For  the  case  that  the  process  noise  Q(k)  =0,  for  all 
k,  we  can  take  4Kk)  to  be  defined  by  • 

WO)  =  P1(0) 

4»(k)  =  Pn  (k)  $T(0,k) 


2] 


5.  DETECTION  AND  ESTIMATION  OF  JUMP  USING  THE  GLR  APPROACH: 
ALGORITHM  I 

From  (41)  we  see  that  Ay  00  takes  on  only  the  values  of  zero 
and  A  y  (q)  where 

A y (q)  -i >~1(q)  Axq  (44) 

It  is  zero  before  the  jump  and  Ay (q)  after  the  jump.  The 
residual  to  be  minimized  after  the  jump  is 

Av  00  =  AZ00  “  AH00  AY  (45) 

which  has  covariance  a R(k)  as  defined  by  (25).  Here,  we  use 
Ay  to  denote  the  unknown  constant.  For  each.k,  we  desire 

A  A  A 

to  obtain  the  estimates  q(k)  and  Ay  (q(k),k)  that  render 
a  minimum  to  the  function 

*3  T  -1 

J(q,  Ay;k)  =  I  lAz(i)l  [tfi(i)l  t a* Ci> ] 
i=l 

k  T  -1 

+  l  lAz(i)  -  AH(i)  Ay]  lAR(i)l  lAz(i)  -  AH(i)  Ay]  (46) 

i=q+l 

or,  equivalently,  a  maximum  to  the  function 
k  T  -1 

«fq»  AyyR)  -  i  [  a z (i) 3  (aRU) ]  tAZ(i)]  -  J(q»  Ay»'k)  (47> 

i=l 

This  latter  fucntion  is  the  logarithm  of  the  generalized 
likelihood  ratio  (GLR),  [5].  Note  that  the  argument  q  appears 
only  in  the  limits  of  the  sum  in  (46) .  In  the  approach  of 
[5]  the  jump  time  appears  in  the  terms  of  the  sum  as  well 
as  in  the  limits  of  the  sum.  For  each  possible  jump  time 

a 

q  the  optimum  value  Ay  (<3/fc)  satisfies 
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(48) 


C(q;k)  Ay(q,k)  =  d(q;k) 
where 

^  T  -1 

C (q;k)  =  I  [A  H  (i) )  [ AR(i) ]  ^ [AH (i) ] 

i=q+l 


d(q;k)  =  £  [AH(i)]T[AR(i)3_1[Az(i)] 

i=q+l 


(49) 

(50) 


In  view  of  the  solution  (48)  ,  the  log  likelihood  ratio  (47) 
simplifies  to 


A(q,  Ay(q,k);k)  =  AyT(q,k)  C(q;k)  Ay(q*k)  *  dT(q;k)^y (q,k)  (51) 

or,  equivalently, 

Jl(q,  Ay(q»k);k)  »  dT(q;k)  cT^q.-k)  d(q;k)  (52) 


The  maximum  likelihood  estimate  of  the  jump  time  is 

A  A 

q (k)  =  arg  max  Mq,  Ay(q,k);k) 

q 

and  the  maximum  likelihood  estimate  of  the  jump  magnitude 


Ay  (q  (k)  ,k)  ■  C_1  (q(k);k)  d(q(k);k) 


(53) 


(54) 


A  jump  is  detected  at  time  k  if  a  threshold  is  exceeded, 

[5]: 

A(q,  Ay(q(k),k))  >  2  in  (n)  (55) 

where  the  value  r\  is  chosen  to  provide  a  reasonable  tradeoff 
between  false  and  missed  alarms. 
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A  difference  between  our  approach  and  that  of  [5]  is 
that  the  failure  signature  matrix  G(j;e)  depends  on  the  jump 
time  and  our  corresponding  matrix  ^H(j)  does  not.*  As  a  result 
of  this  advantage  the  computation  of  updating  several  sequences 
of  matrices  is  avoided.  For  the  purpose  of  estimating  the 
jump  time  and  the  jump  magnitude  at  time  k  it  suffices  to 
have  the  following  matrices  in  storage: 

C  (  P-  i)  ,  i  =  1,2  ....  k 
d (  0; i)  ,  i  =  1,2,  . . .  fk 

The  matrices  C(0;k)  and  d(Q;k)  are  computed  recursively 

C (  0; k)  =  C(0;k-1)  +  /£T(k)  tfl(k)”1  &( k)  (56) 

d  (  0;  k)  =  d  (0i;k-l)  +  AHT(k)  ARW"1  A  z  (k)  (57) 

The  matrices  C(q;k)  and  d(q;k)  are  obtained  by  subtraction 

C (q;k)  =  C  (  0; k)  -C(0;q),  0<  q  <k  (58) 

d(q;k)  =  d(0?k)  -  d  (Q?q),  0<q<k  (59) 

It  is  unnecessary  to  evaluate  (52)  for  all  q,  0<q  <k.  A 
simple  search  procedure  can  be  employed  to  locate  the  maximizing 

A  a  A 

q(k)  of  (53):  After  obtaining  A  y  we  use  (40)  to  estimate 

A*(k)  (k)  Ay  (60) 

♦This  is  discussed  in  the  next  section. 

*In  the  search  procedure  one  may^use  the  Gaussian  elimination 
method  to  compute  the  solution  Ay(<3*k)  of  (48)  at  some  o<q<k. 
Therefore,  the  log  likelihood  ratio  (51)  may  be  evaluated 
without  inverting  £(q;k).  This  matrix  need  only  be  inverted 
at  the  maximizing  <^(ki,  provided  (55)  is  satisfied,  to  obtain 
the  covariance  of  Ay(q(k),k). 
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Its  covariance  is  given  by 

AP(k)  =  E{A  e (k) A  e(k)T}  -  ^ (k)  C_1  (q (k) ;k)^T  (k)  (61) 

A 

where  Ae(k)  =  Ax(k)  -  Ax(k).  The  covariance  of  the  estimate 

A 

Ay  is  given  by 

E{  A  e1(k)  Ae1(k)T}  =  C_1  (q(k);k)  (62) 

A  A 

where  A©! (k)  =  ayOO  "  AY  (q(k) ,k) .  This  covariance  is  based 

on  the  assmption  that  the  a  priori  covariance  for  aY  time 

qlk)  is  infinite,  [16,  p.  206].  The  estimates  given  by  (60) 

and  (61)  are  used  as  starting  values  in  the  A~filter.  Eqs.  (28)  and 

(29)  are  used  to  reinitialize  the  filter  of  for  the  case  of  multiple  jumps. 

If  the  a  priori  covariance  for  Ax(q)  is  not  infinite 
but  is  given  by  AP (q)  with  inverse  AP-1  (q)  the  a  priori 
covariance  AP.^  (q)  for  Ay  (q)  has  inverse 

AP1”1(q)  =  ipT(q)  AP“1(q)  i|>(q)  (63) 

A 

If  Ax(q)  is  the  mean  of  Ax (q) , the  mean  of  Ay(q)  is 
Ay  (q)  =  i|>  1  (q)  Ax(q) 

With  the  modifications 

C*  (q;k)  =  C  (q;k)  +  APj^”1  (q) 
d*  (q;k)  *  d(q;k)  +  A?!*1 (q)  AY<<l) 

The  method  of  solution  is  as  outlined  above. 


(64) 

(65) 

(66) 
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If  we  know  that  the  jump  is  caused  by  a  lower  dimensional 


vector  such  as  &uq  where 

AXq  =  D  4uq  (67) 

we  proceed  as  follows.  From  (44)  and  (67)  we  have 

Ay(q)  =^-1(q)  DAuq  (68) 

we  define 

C(q;k)  =  D%  ~T(q)  C  (q;k)  ty  _1  (q)D  (69) 

d (q;k)  =  D%"T(q)  d(q;k)  (70) 


and  use  the  method  of  solution  as  described  above. 

The  GLR  technique  described  in  [9]  requires  the  implementation 
of  a  Kalman-Bucy  a -filter  (employing  the  computational  savings 
techniques  discussed  therein  for  each  observation  time  in 
the  past.  Several  matrices  have  to  be  stored  and  updated 
for  each  observation  time  in  the  past.  The  attractive  feature 
of  that  technique  is  that  it  requires  neither  the  inverse 
of  C(q;k).nor  the'solution  to  (48).  Instead  it  requires  at  most 
the  inverse  of  a  matrix  having  the  dimensions  of  the  largest  .  ' 
correlation  block  of  the  measurement  noise  covariance  matrix; 
sequential  updating  of  the  components  .of  the  measurement  vector  ‘ 
is  used.  No  matrix  inversion  is  required  by  the  Kalman-Bucy  filter 
for  uncorrelated  measurement  noise.  We  use  the  sequentially  updated 
Kalman  filtering  technique  together  with  the  decomposition  of  the 
maneuver  signature  matrix  to  develop  our  algorithm  II  in  Section  7. 
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6.  COMPARISON  WITH  PREVIOUS  GLR  TECHNIQUES 

The  a  posteriori  jump  formulation  is  used  in  (4)  and  (5). 

It  has  the  following  relationship  with  the  a  priori  jump  formulation. 
Define : 


p  =  q+i 

(71) 

AXp  =  4>(q+l,q)  AXq 

(72) 

Equation  (4)  can  be  rewritten  as: 

x  (k+1)  *  d>(k+l,k)  x  (k)  +  Axp  5pfk+i  +  r(k)  w(k)  (73) 

where  p  is  the  jump  time  and  Axp  is  the  jump  magnitude.  The 
effect  of  the  jump  appears  in  the  state  x(p)  and  in  the  measurement 
z(p)  at  the  jump  time;  that  is, it  occurs  right  before  the  measurement 
z(p).  This  formulation  is  used  in  [3]  -  [9],  Therein,  the 
effect  of  the  jump  on  the  innovations  is  analyzed.  The  general 
form  is  given  by,  [3]  and  [5], 

Yl(k)  =  G(k;p)  Axp  +  Yi(k)  (74) 

where  Yi(k)  and  Y^(k)  are  the  predicted  measurement  residuals 
under  Conditions  Hi  and  Hi,  respectively.  The  matrix  G(k;p) 
is  called  the  failure  signature  matrix+  and  it  is  computed  by 
the  following  recursive  algorithm,  [5]  and  [6].  At  each  observation 
time  k  the  following  matrices,  having  been  computed  previously 
at  time  k-1, are  held  in  storage: 

■•■Herein,  we  refer  to  it  as  the  maneuver  signature  matrix. 
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G(k-l;j),  k-l-M  <j  <k-l 


(75) 


F(k-l;j),  k-l-M  <j  <k-l  (76) 

where  it  is  assumed  that  a  sliding  window  of  length  M  is  being 
used  to  detect  and  estimate  the  jump.  For  each  k  the  matrix 
G(k;k)  satisfies 

G(k;k)  =  H (k) ,  all  k  (77) 

The  following  matrices  are  computed  at  the  time  k: 


$(k; j) 

*  $  (k ,  k 

-1)  *(k 

-l,j),  k-M<j<k 

(78) 

G(k;j) 

-  H (k ) [*(k,j)  - 

S  (k; j ) 1 ,  k-M<j<k 

(79) 

S (k  ; j ) 

*  4>(k,k- 

-1) F (k-1 

« j ) *  k-M<j<k 

(80) 

F  (k  ;  k ) 

=  KX  (k) 

H  (k) 

(81) 

F (k ; j  ) 

-  KX(k) 

G(k;j) 

+  S  (k; j ) ,  k-M£j<k 

(82) 

The  following  matrices  have  also  been  computed  and  stored  at 
the  previous  observation  time  k-1: 


^w 

(k-l;j),  k-Mcj^k-l 

(83) 

<*w 

(k-l;j),  k-M<j<k-l 

(84) 
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At  the  new  observation  time  k  they  are  updated  using  the  equations 


Cw  (k;k)  *  Ht (k)  Vi-1<k)  H (k)  (85) 

dw  (k ;  k )  =  HT(k)  Vx-^k)  Yx(k)  (86) 

Cw  ( k ;  j )  =  GT(k;j)  Vi'  ;k)  G(k;j)  +  Cw(k-l;j),  k-M<j<k  (87) 

dw(k;j)  *  GT(k;j)  Vi”1^)  ^(k)  +dw(k-l;j),  k-M<j<k  (88) 

The  log-likelihood  ratio  to  be  maximized,  [5]  and  [6] ,  is 

J,w(k;  j)  =  dwT  (k ;  j )  Cw-l(k;j)  dw(k;j),  k-M<j<k  (89) 

which  is  evaluated  at  past  observation  times  j  in  order  to  determine 
its  maximum  value  and  compare  it  to  a  set  threshold  for  jump 
detection. 


In  our  approach  Eq.  (74)  is  given  by 

Yl(k)  -  H (k)  $(k,k-l)  AX(k-l)  +  ^(k)  (90) 

since 


Yl(k)  =  z(k)  -  H(k)  <j>(k,k-l)  xi(k-l) 
Y2(k)  =  z (k)  -  H (k)  $(k,k-l)  x2(k-l) 
Yl(k)  -  Y2<*> 


From  Eqs.  (18)  and  (19)  we  have 


Ax(k-l)  ■  A$(k-l.fP)  Ax  (P) 


AX(p)  ■  [I-Ki(p)  H (p) ] AXp 
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(91) 

(92) 

(93) 


(94) 

(95) 


Comparing  (74)  and  (90)  we  see  that 


G(k;p)  Axp  =  H (k )  *(k,k-l)  A*(k-l,p)  [I-Ki(p)H(p) ) *xp  (96) 
Using  (38)  we  rewrite  (96)  as 

G(k;p)  Axp  =  H(k)  4>(k,k-l)  ip(k-l)  +_1(p)  [I-Kx(p)  H(p)]Axp 

(97) 

Consequently,  in  our  approach  G(k;p)  ^xp  is  decomposed  as 

G(k;p)  ^xp  »  Gi(k)  G2(p)  Axp  (98) 

where 

Gi(k)  =  H (k)  $(k,k-l)  ^ (k-1)  (99) 

G2(P)  *  ^_1(P)  [I-Ki(p)  H (p) ]  (100) 

Using  (18),  (71),  and  (72)  we  have 

G2(P)  Axp  *  \Jj-^(g)  Axq  (101) 

Using  (44),  Eq.  (101)  becomes 
G2  (P)  Axp  "  Ay(q) 

Consequently,  Eq.  (98)  can  be  written  as 

G(k;p)  Axp  *  Gi(k)  AY (P-1)  (102) 
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The  decomposition  given  in  Eq.  (102)  results  in  the  matrix 
Gi(k)  which  does  not  depend  on  the  jump  time  and  in  the  vector 
Ay(p-l)  which  satisfies  the  constant  state  equation  (41). 

The  computational  burden  of  the  GLR  technique  of  [5]  and  (6] 
is  directly  related  to  the  dependency  of  G(k;p)  on  the  jump 
time.  Because  of  this  dependency  the  matrix  G(k;j)  must  be 
computed  and  stored  for  each  candidate  jump  time  j  in  the 
past.  This  requires  at  each  new  time  k  the  computation  and 
storage  of  the  matrices  given  in  (75)  -  (88)  for  all  j  in 
a  sliding  window  k-M  <j_<k  of  candidate  jump  times. 

In  our  approach  we  are  actually  using  ^l(k)  as  given  in 
(42)  rather  than  G]_(k)  as  given  by  (99).  This  is  because 
we  are  using  the  a  posteriori  measurement  residuals  ^z(k) 
given  in  (20)  rather  than  the  a  priori  measurement  residuals 
Yl(k)  given  in  (91).  That  is,  ^H(k)  is  the  resulting  decomposition 
matrix  for  the  a  posteriori  measurement  residual  approach. 

An  analysis  of  the  GLR  algorithm  requirements  for  storage 
and  multiplications  per  storage  update  is  given  in  Table  5 
for  the  approach  of  [5]  and  [6]  and  in  Table  6  for  our  approach.* 
Let  the  storage  requirements  be  denoted  by  Msw  for  that  of 
[5]  and  [6]  and  by  Ms  for  our  approach.  The  storage  requirements 
are 


*The  measurement  matrix  R(k)  is  assumed  to  be  diagonal  and 
all  elements  of  the  matrices  H(k)  and  $(k,k-l)  are  multiplied 
in  matrix  products. 
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TABLE  5 


WILLSKY  AND  JONES'  GLR  ALGORITHM  REQUIREMENTS: 
STORAGE  AND  MULTIPLICATIONS  PER  STORAGE  UPDATE 

k  =  CURRENT  OBSERVATION  TIME 
n  =  DIMENSION  OF  STATE  VECTOR  x 
m  =  DIMENSION  OF  MEASUREMENT  VECTOR  z 
M  =  LENGTH  OF  SLIDING  WINDOW 


I .  STORAGE  REQUIREMENTS 


MATRICES 

DIMENSION 

STORAGE  RE< 

1. 

{G(k; j) 

:k-M< j<k} 

mxn 

(M-l) nm 

2. 

{  F  ( k ;  j ) 

:.k-M<  j  <k  } 

nxn 

Mn2 

3. 

{Cw(k; j 

) :k-M<j<k} 

nx (n+1) 

2 

Mn (n+1) 
2 

4. 

{dw ( k ; j 

)  : k-M< j  <k } 

nxl 

Mn 

5. 

{4>(k;  j) 

:k-M<j<k} 

nxn 

(M-l) n2 

TOTAL  STORAGE  =  (M- 

-l)Miw  +  M2w 

wher^ 

u  r5n2  3n  , 

Mlw  *  [— +  nm  +y~  ] 


«2w  -  + 


1“  1 
2  J 


Table  5  (continued) 


II.  MULTIPLICATIONS  NEEDED  TO  UPDATE  STORED  MATRICES 


EQUATIONS 

NUMBER  OF 
MULTIPLICATIONS 

1. 

Eqs . 

(79) 

and 

(80) 

for  G(k;j),  k-M<j<k 

(M-l) [n3  +  n2m] 

2. 

Eqs . 

(81) 

and 

(82) 

for  F(k;j);  k-nuj^k 

M[n2m] 

3. 

Eqs . 

(85) 

and 

(87) 

for  Cw(k,j),  k-m<j£k 

M  [  2  +  nmz] 

4. 

Eqs . 

(86) 

and 

(88) 

for  dw(k, j) ,  k-M< j <k 

M[nm  +  m2] 

5. 

Eq. 

(78) 

for 

$(k»  j) 

,  k-M<j<k 

(M-l)  [n3] 

TOTAL  MULTIPLICATIONS  =  (M-2)  Niw  +  N2W  where 

N]_w  =  [2n3  +  n2m  +  nm2  +  nm  +  m2] 
n2w  =  l2n3  +  4n2m  +  2nm2  +3nm  +  2m2^ 


TABLE  6 


NEW  GLR  ALGORITHM  I  REQUIREMENTS 
STORAGE  AND  MULTIPLICATIONS  PER  STORAGE  UPDATE 

k  =  CURRENT  OBSERVATION  TIME 
n  =  DIMENSION  OF  STATE  VECTOR  x 
m  =  DIMENSION  OF  MEASUREMENT  VECTOR  z 
M  =  LENGTH  OF  SLIDING  WINDOW 


I.  STORAGE  REQUIREMENTS 


1. 

2. 

3. 

4. 

5. 

6. 


MATRICES 

DIMENSION 

STORAGE  REQUIREMENTS 

\p  ( k ) 

nxn 

n2 

A«  (k ) 

mxn 

nm 

AZ  (k) 

mxl 

m 

AR(k)-l 

mxm 

m (m+1) 

2 

{C  ( 0 ;  j  )  :k-M<j  <k} 

nxn 

..n  (n+1) 

M  2 

{d  (0;  j  )  :k-M<j<k} 

TOTAL  STORAGE  =  (M-l) 
„  r  02  ,  3n  , 


Mo  -  t3n2 

"2  "  l~ 


+  nm  .  + 


nxl 

Mi  +  M2  where 

n»2  3n  ,  3ro  . 

2  2~  2  i 


Mn 
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Table  6  (continued) 


II.  MULTIPLICATIONS  NEEDED  TO  UPDATE  STORED  MATRICES 


NUMBER  OF 


EQUATIONS 

MULTIPLICATIONS’ 

1.  Eq. 

(38) 

for 

i|)(k) 

n2 

2.  Eq. 

(42) 

for 

AH(k) 

n2m 

3 .  Eq . 

(20) 

for 

Az  (k) 

nm 

4.  Eq. 

(25) 

for 

AR(k)-1 

f  "2 

+ ! 

5 .  Eq. 

(56) 

for 

C  ( 0  ;  k  ) 

n2m 

2 

tnm2.SE 

6 .  Eq . 

(57) 

for 

d(0;k) 

nm  + 

m2 

TOTAL  MULTIPLICATIONS  *  (M-2)Ni  +  N2  where 


Ni  -  0 

„  r  t  .  3n2m  .  _ o  ,  5nm  ,  5m2  ,  m  , 

N2  =  [nJ  +  — ^ —  +  nm*  +  -y  +  — —  +  j  J 


*For  the  case  that  the  process  noise  is  zero  it  is  not  necessary 
to  compute  \p(k)  unless  a  jump  isTdetected.  It  suffices  to  compute 
AH (k)  where  AH (k)  =  H(k)  P. (k)  $  (0,k).  This  product  is  performed 
with  at  most  2n2m  multiplications  which  replaces  the  sum  n3  +  2n'Tn  . 
If '  A '  jump-" is-  detected  we.  need  \p tk)  for  computing  A*(k)  from 
Ay  (k) . 


(M-l)  M^w  +  M2w 


M  sw  * 


where 

Mlw 

**2w 
and  are 


5n2 


+  nm  + 


3n 

T 


3n2  3n 

2  2 


Ms  =  (M-l)  Mx  +  M2 


where 

Mi 


n2  3n 
2  +  2 


Mo  = 


3n2  +  nm  4-  m2  4.  3n.  3m 
—  +  nm  +  -j-  +  j-+  y 


Consider  the  differences 


AM^  =  M^w  -  Mj_  =  2n^  + 


nm 


a  _  m2  3m 

AM2  =  M2w  -  M2  =  -nm  -  —  ~  2~ 


Define 


AMg  *  Mgw  -  Mg 


(103) 

(104) 

(105) 

(106) 

(107) 


(108) 

(109) 


(110) 


as  the  difference  in  storage  requirements.  The  dimension 
m  of  the  measurement  vector  is  usually  much  less  than  the 
dimension  n  of  the  state  vector.  For  sliding  windows  with 
_ M  >1  we  have  the  inequality 
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A  Mg  >(M-2)  (n2  +  nin)  (HD 

Let  the  number  of  multiplications  be  denoted  by  Nw  for 
the  approach  of  [5]  and  [6]  and  by  N  for  our  approach. 

The  number  of  multiplications  needed  to  update  the  stored 
matrices  are 

Nw  *  (M-2)  Nlw  +  N2w  (112) 

where 

N]_w  =  2n3  +  —  ~~  +  nm2  +  nm  +  m2 

N2w  =  2n3  +  4n2m  +  2nm2  +  3nm  +  2m2 

and 

N  =  (M-2)  Ni  +  (115) 

where 

NL  =  0  (116) 

m  ^  ,  3n 2m  j  ....  9  i  5nm  ,  5m2  ,  ro  ^  ^  i  ^  \ 

N2  =  nJ  +  — ^ —  +  nm*  +  — —  +  — j-  +  (117) 

Eq.  (116)  shows  that  our  GLR  algorithm  requires  no  multiplications 

to  update  stored  matrices  at  past  observation  times  j, 

j<k. 


(113) 

(114) 
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Consider  the  differences 


I 


p  ! 


i 


ANi  =  Niw  -  Ni  =  Niw 


(118) 


AN2  =  N2w  “  ^2 


=  n3 


5n2m 


+  nm 


2  + 


nm 


m2 


m 

2 


(119) 


Define 

AN  =  Nw  -  N  (120) 

as  the  difference  in  the  number  of  multiplications. 

For  sliding  windows  with  M  >2  we  have  the  inequality 

AN  > (M-2)  (n3  +  2n2m)  (121) 

The  inequality  (111)  demonstrates  the  savings  in 
storage  provided  by  our  GLR  technique.  The  inequality 
(121)  demonstrates  the  savings  in  multiplications.  These 
savings  are  a  direct  result  of  the  decomposition  discussed 
above  and  given  in  Eq.  (102)  for  the  a  priori  measurement 
residual  approach  and  given  by 

H(k)  ^(k)  Ay(q)  (122) 

for  the  a  posteriori  measurement  residual  approach. 
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Both  of  the  above  GLR  techniques  require  the  inverse-4-  of 
an  nxn  matrix  in  order  to  evaluate  the  log-likelihood  ratio 
(52)  or  (89) .  Since  our  approach  requires  no  updating 
of  matrices  at  past  observation  times  it  suggests  using 
an  innovations  based  scheme  to  determine  the  current  observation 
times  k  at  which  one  should  look  for  a  jump  in  the  past. 

If  the  innovations  appear  to  be  "normal"  no  jump  in  the 
past  is  to  be  searched  for  by  evaluating  (52);  consequently, 
no  inverse  is  taken.  It  is  when  the  innovations  appear 
to  be  less  than  "normal"  that  a  search  is  to  be  made  for 

A 

the  optimizing  jump  time  q  of  (52).  We  now  discuss  the 
GLR  technique  of  [8]  and  [9]  which  requires  no  inverse 
of  an  nxn  matrix.  We  treat  the  impulse  input  version  of  [8]  and 
[9]  rather  than  the  step  input  case;  see  Appendix  E. 

The  GLR  technique  [9]  requires  a  Kalman-Bucy  A-filter 
for  each  observation  time  j,  k-M<j£k  where  k  is  the  current 
time  and  M  is  the  length  of  the  sliding  window.  In  our  notation 
and  for  the  jump  system  (73)  and  (5)  that  algorithm  utilizes 
the  maneuver  signature  matrix  G(k,j)  in  the  form 

G(k,j)  =  H (k)  $(k,k-l)  Aj(k-l),  k-M<j <k  (123) 

where  Aj(k)  is  given  by 

Aj(j)  =  [I  -  Ki(j)  H ( j ) ]  (124) 

“  A 

It  suffices  to  solve  (48)  for  Ay(<3rk).  The  inverse  C_1(q;k) 
needs  only  to  be  computed  at  the  maximizing  argument  of  (53) 
when  (55)  is  satisfied.  Since  jumps  are  infrequent  (55)  will 
be  satisfied  only  infrequently.  Consequently,  inverses 
are  seldom  required. 
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I 


Aj (k)  =  A*  (k , k-1)  Aj (k-1) 


(125) 


We  make  the  definitions 

AHk(k)  =  H(k)  (126) 

AHj(k)  =  H(k)  4>  (k,k-l)  Aj  (k-1) ,  j<k  (127) 

At  each  observation  time  k  the  following  matrices,  having 
been  computed  previously  at  time  k-1,  are  held  in  storage: 


Aj  (k-1) , 

k-l-M< j<k-l 

(128) 

A 

Avj (k-1) , 

,  k-l-M<  j<.k-l 

(129) 

APj(k-l), 

,  k-l-M<j<k-l 

(130) 

dj (k-1) , 

k-l-M<j<k-l 

(131) 

where  it  is  assumed  that  a  sliding  window  of  length  M  is  being 
used  to  detect  and  estimated  the  jump.  At  the  current  time  k 


the  following  matrices  are  computed: 

Ak(k)  =  [I  -  Ki(k)  H (k) 3  (132) 
Aj  (k)  =  A4>(k,k-l)  Aj  (k-1) ,  k-M<j<k  (133) 
AHj(k)  *  H(k)  $(k,k-l)  Aj  (k-1 )  ,  k-M<  j<k  (134) 
AVj(k)  *  AHj (k)  APj(k-l)  AHjT(k)  +Vi(k),  k-M< j<k  (135) 
AKj(k)  *  APj(k-l)  AHjT(k)  AVj_1(k),  k-M<j<k  (136) 
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9 


AVj(k)  *  Avj  (k-1)  +  AKj  (k)  fYl(k)  -  AHj  (k)  Av  j  (k-1)  ] 


k-M<j£k  (137) 

APj(k)  =  [I  -AKj(k)  AHj (k)  ]  APj  (k-1)  ,  k-M<j<k  (138) 

dj(k)  -  dj(k-l)  +  AHjT(k)  V!"1  (k)  Yi(k),  k-M<j<k  (139) 

4  j  (k)  =djT(k)  AVj  (k)  ,  k-M<j  <k  (140) 

where 

d^  (k-1)  *  0  all  k  (141) 

AVj^tk-l)  ■  0  all  k  (142) 

APk“1(k-l)  =0  all  k  (143)+ 


The  above  equations  constitute  a  Kalman-Bucy  a -filter 
for  each  candidate  jump  time  j  in  the  window  k-M<j<k.  The 

/v 

estimate  AVj (k)  is  the  optimal  estimate  at  time  k  of  Axp  when 
the  jump  time  p  coincides  with  the  candidate  jump  time  j. 

A 

If  the  jump  threshold  is  exceeded  at  the  maximizing  j  *  p 
of  (140)  the  optimal  estimate  of  x(k)  is  given  by,  [9], 

x (k)  =  xx(k)  +  Ap (k)  Avp ( k )  (144) 

and  its  covariance  by 

P  (k)  =  Pi(k)  +  Ap (k)  A  Pp (k)  A^(k)  (145) 


+The  initial  covariance  APk(k-l)  is  defined  as  the  identity 
matrix  I  times  a  very  large  number. 


41 


The  sequentially  updated  Kalman  filtering  technique,  [10] 
and  [llj,  is  used  in  [8]  and  [9]  to  avoid  the  matrix  inverses 
in  (136)  and  (139)  for  the  case  when  some  of  or  all  the  components  of 
the  measurement  vector  are  uncorrelated.  Note  that  the  maximizing 
jump  time  j  is  easily  obtained  from  (140),  requiring  only 
Mn  multiplications. 


An  analysis  of  the  GLR  algorithm  requirements  for  storage 
and  multiplications  per  storage  update  is  given  in  Table  7 
for  the  approach  of  [8]  and  [9] .  Let  the  storage  requirements 
be  denoted  by  Msc.  The  storage  requirements  are 

Msc  =  (M-l)Mic  +  M2c  ( 146 ) 


where 


(147) 


+  2nm  +  m2  +  — 


(148) 


The  differences  in  storage  requirements  between  the  above 
approach  and  our  approach  are: 

AMic  =  Mic  --Mi  =  n2  +  n  (149) 

m2  3m 

AM2c  =  M2c  -  Ml  =  nm  +  y-  +n-jL  (150) 
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Define 


AMgc  =  Mgc  —  Mg  (151) 

as  the  difference  in  storage  requirements  for  the  two  approaches. 
For  sliding  windows  with  M>1  we  have  the  inequality 

AMsc  > (M-l)  [n2  +  n]  (152) 

Let  the  number  of  multiplications  be  denoted  by  Nc  for 
the  approach  of  [9] .  The  number  of  multiplications  needed 
to  update  the  stored  matrices  are 

Nc  =  (M-2)  Nlc  +  N2c  (153) 

Where 

Nic  =  +  5  n2ro  +  n2  +  4nm  +  n+m  (154) 

N2c  =  n2  +  9n2m  +  2n2  +  8nm  +  2n  +  2m  (155) 

Define 

AN2c  =  N2c  -  N2  (156) 

We  have  the  inequality 

&N2c  >  6  n2ro  +  2n2  +  3nm  (157) 
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It  is  not  fair  to  make  a  direct  comparison  between 
N]_c  and  Ni  since  the  approach  of  [9]  provides  the  estimate 

A 

Avj(k)  for  each  j  in  the  sliding  window  and  ours  does  not. 

In  order  to  make  a  fair  comparison  we  must  add  in  the  number 

/* 

of  multiplications  needed  to  solve  (48)  for  Ay(q»k).  Since 
C(q;k)  is  a  symmetric  matrix,  the  Cholesky  method,  [18], 
may  be  employed  to  solve  (48)  ;  the  number  of  multiplications 


required  are 


N,  _  n 2  .  3n£  n 

**3  -  g-  +  +  3 


(158) 


Since  we  employ  a  search  procedure  to  maximize  (53) 
it  is  not  necessary  to  solve  (48)  at  each  q  in  the  sliding 
window.  Let  us  assume  at  the  very  worst  that  we  will  need 
to  evaluate  (51)  at  M-2  observation  times  in  the  sliding 
window;  that  is,  we  need  not  solve  (48)  at  two  points.  Define 


ANlc  =  Nic  -  N3 


We  have  the  inequality 


ANic  n^  +  4n^ 


Define 


A  Nc  =  (M-2) 4  Nlc  +  4  N2c 


(159) 


(160) 


(161) 
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For  sliding  windows  with  M>2  we  have  the  inequality 

ANc  >  (M-2)  +  4n2]  (162) 

The  inequalities  (152)  and  (162)  demonstrate  the  savings 
in  storage  and  in  multiplications  provided  by  our  GLR  technique 
as  compared  to  that  in  [9] . 


TABLE  7 


CHANG  AND  DUNN'S  GLR  ALGORITHM  REQUIREMENTS:* 
STORAGE  AND  MULTIPLICATIONS  PER  STORAGE  UPDATE 

k  =  CURRENT  OBSERVATION  TIME 
n  =  DIMENSION  OF  STATE  VECTOR  x 
m  =  DIMENSION  OF  MEASURMENT  VECTOR  z 


M  =  LENGTH 

OF  SLIDING 

WINDOW 

I. 

STORAGE  REQUIREMENTS 

MATRICES 

DIMENSION 

STORAGE  REQUIREMENTS 

1 

.  {Aj  (k)  :  k-M<j£k  } 

nxn 

Mn2 

2 

.  (Avj(k):  k-M<j<k} 

nxl 

Mn 

3 

.  {  APj  (k)  :  k-M<j  <k  } 

nxn 

wn (n+1) 

M  2 

4 

.  {dj  (k)  :  k-M<j  <k  } 

nxl 

Mn 

5 

.  AHj(k)  -  DUMMY  MATRIX 

mxn 

nm 

6 

.  AVj(k)  -  DUMMY  MATRIX 

mxm 

m2 

7 

.  AKj(k)  -  DUMMY  MATRIX 

nxm 

nm 

TOTAL  STORAGE  =  (M-l) 

Mic  +  m2c  where 

„  3n2  5n 

Mlc  =  2  +  2 

M  3n2  5 

M2c  ~  — 2~  +  2nm  +  m^  + 

5n 

2 

♦Impulse  input  case. 
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TABLE  7  (continued) 


II.  MULTIPLICATIONS  TO  UPDATE  STORED  MATRICES 

EQUATIONS  NUMBER  OF  MULTIPLICATIONS* 


1.  Eqs. 

(132) 

and  (133) 

for  Aj  (k)  , 

(M-l)  n3  + 

k-M<j  <k 

2 .  Eq . 

(137) 

A 

for  A  v j  { k )  , 

k-M<j <k 

M2nm 

3.  Eq. 

(138) 

for^P j  (k)  , 

k-M<  j<k 

M2n2m 

4.  Eq. 

(139) 

for  dj (k)  , 

k-M< j£k 

M[nm  +  m] 

5 .  Eq. 

(134) 

for  a  Hj  (k)  , 

k-M< j< k 

(M-l) 2n2m 

6.  Eq. 

(135) 

for  AVj (k)  , 

k-M<j<k 

M[n2m  +  nm] 

7.  Eq. 

(136) 

for  AKj  (k)  # 

k-M<  j<k 

M [n2  +  n) 

TOTAL  MULTIPLICATIONS 

=  (M-2)  Nlc  + 

N2c  where 

Nlc 

=  n3  + 

5n2m  +  n2 

+  4nm  +  n+m 

N2c 

=  n3  + 

9n2m  +  2n2 

+  8nm  +  2n  + 

2m 

*The  sequentially  updated  Kalman  filtering  technique,  [10]  and 
[11] ,  is  used. 


7. 


BANK  OF  A-FILTERS  USING  DECOMPOSITION  OF  THE  MANEUVER 


SIGNATURE  MATRIX:  ALGORITHM  II 


The  GLR  technique  [9]  which  uses  the  maneuver  signature 
matrix  defined  by  (123)  -  (125)  requires  the  computation  and 
storage  of  the  matrix  Aj(k)  for  each  observation  time  j  in 
the  past.  From  Table  7  we  see  that  this  requires  the  following 
number  of  multiplications: 

N  =  (M-l)  n3  +  n2m 

cl 

for  sliding  windows  of  length  M.  From  Table  7  we  note  that 
this  matrix  is  the  only  one  which  requires  multiplications 
proportional  to  the  third  power  of  n.  It  follows  that  N 

cl 

is  proportional  to  n4  when  M>n. 

The  dependence  of  Nfl  on  M  can  be  removed  by  using  our 
decomposition  (42)  or  (122)  of  the  maneuver  signature  matrix. 

We  develop  this  idea  next. 

Consider  the  employment  of  a  bank  of  Kalman-Bucy  constant 
A-state  filters  for  a  moving  window  of  length  M.  That  is, 
at  each  j,  k-M<j<k,  we  employ  the  constant  A-state  filter 
defined  by  (41)  -  (43)  .  The  filtering  equations  are  given 
in  Table  8.  The  GLR  algorithm  is  as  follows. 


(163) 
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TABLE  8 

FILTERING  EQUATIONS  FOR  THE  CONSTANT  A -STATE  SYSTEM 

(j  =  CANDIDATE  JUMP  TIME) 


A  Yj  (k+1 1  k)  =  AYj  (k)  »  k>j-l  (i) 

A  P j  (k+1 1  k)  =AP j  (k)  ,  k>j-l  (ii) 

AYj  (k)  =  Az<k)  -  AH(k)  A Y j  0^ | )c-l )  ,  k >j  (iii) 

A  Vj  (k)  =  A  H  (k)  APj  (k  |k-l)  AHT(k)  +  AR  (k)  ,  tej  (iv) 

AKjtk)  =  A  Pj  (k|k-l)AHT(k)  AVj"1(k),  k>p  (v) 

/Ty  •  (k)  -  Ay.  (k|k-l)  +AK.|k)  AY.  (k)  ,  (vi) 

J  J  J  J 

A  P  j  (k)  =  [I  -  AKj(k)  A  H  (k)  ]  APj(klk-l),  k>j  (vii) 

where  k  =  current  observation  time  and 

A  H (k)  =  H(k)  Y  (k) ,  all  k  (viii) 

A  R(k)  -  R(k)  V1_1(k)  R(k) ,  all  k  (ix) 

A  z  (k)  =  z(k)  -  H(k)  x^k),  all  k  (x) 

APj(j-l)  =  c  I,  c  a  very  large  number  (xi) 

AYj(j-l)  =  0  (xii) 

y  (0)  =  I  (xiii) 

y  (k)  *  a*  (krk-l)  y  (k-1) ,  k>  0  (xiv) 

A*  (k,k-l)  «  [I  -  K^k)  H(k)  ]  #(k,k-l)  ,  k>  0  (xv) 
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At  each  observation  time  k  the  following  matrices,  having 


been  computed  previously  at  time  k-1, are  held  in  storage: 

<Hk-l)  (164) 

A 

Ayj(k-l),  k-l-M< j<k-l  (165) 

APj(k-l),  k-l-M<j<k-l  (166) 

dj  (k-1) ,  k-l-M<j<k-l  (167) 

At  the  current  time  k  the  following  matrices  are  computed: 

y  (k)  =  *  (k,k-l)  (k-1)  (168) 

A  H  (k )  =  H(k)  'i'(k)  (169) 

AVj(k)  *  A  H  (k)  APj(klk-l)  AHT(k)  +  AR(k),  k-M<j<k  (170) 

AKj(k)  =  A  Pj  (k  |k-l)  AHT(k)  A  Vj-l(k)  ,  k-M<j<k  (171) 

Ayj(k)  =  A^yj  (k  1  k-1)  +AKj(k)  k-M<j<k  (172) 

A  Pj  (k)  -  [I  -  a  Kj  (k)  A  H(k)  ]  APj(k-l),  k-M<j<k  (173) 

dj  (k)  =  dj  (k-1)  +  AHT(k)  AR_1(k)  Az  (k),  k-M<j<k  (174) 

A 

fc  j(k)  =  djT(k)  Ayj  (k)  ,  k-M<j<k  (175) 

Ak 

A  yk  (k-1)  *  0  (176) 

A (k-1)  *  c  I,  c  a  very  large  number  (177a) 

(177b) 


djt  (k-1)  =  0 


[ 


The  estimate  Ayj(k)  is  the  optimal  estimate  at  time  k 
of ¥_1 (q) Axq  when  the  candidate  jump  time  j  coincides  with 
the  observation  time  q+1.  Note  that  the  filter  for  Ayj  is 
initiated  by  assuming  that  the  a  posteriori  jump  Axq  occurs 

/v 

at  time  j-1.  If  the  maximizing  j  of  (175)  satisfies  (55) 

A 

the  optimal  estimate  Ax(k)  is  computed  as 
Ax(k)  =  'Ft  k )  A  yj  ( k ) 


and  its  covariance  as 

AP(k)  «  ¥  (k)  AP ^ (k)  ^T(k) 

An  analysis  of  the  above  GLR  algorithm  requirements  for 
storage  and  multiplications  per  storage  update  is  given  in 
Table  9.  The  analysis  assumes  the  sequentially  updated  Kalman 
filtering  technique,  (10)  and  [11] ,  is  used.  Let  the  storage 
requirements  be  denoted  by  Md.  The  storage  requirements  are 

Md  =  (M-l)  Mld  +  M2d 


where 


+ 


5n 

2 


n2 

M2d  *  T~ 


+  2nm  +  m2  +  -r 


5  n 


(178) 

(179) 


(180) 

(181) 

(182) 
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TABLE  9 


REQUIREMENTS  OF  NEW  GLR  ALGORITHM  II  -  BANK  OF  CONSTANT  A- STATE 
FILTERS  USING  DECOMPOSITION  OF  MANEUVER  SIGNATURE  MATRIX: 
STORAGE  AND  MULTIPLICATIONS  PER  STORAGE  UPDATE 

k  *  CURRENT  OBSERVATION  TIME 
n  =  DIMENSION  OF  STATE  VECTOR  X 
m  *  DIMENSION  OF  MEASUREMENT  VECTOR  Z 
M  *  LENGTH  OF  SLIDING  WINDOW 


C .  STORAGE  REQUIREMENTS 

MATRICES 

1.  y(k) 

2-  UYj<k)  tk-M<j<k) 

3.  {APj(k):  k-M< j<k} 

4.  {d.(k):  k-M<j<k} 

5.  AH(k) 

6.  AVj Ck) 

7.  AKjtk) 


DIMENSION 

nxn 

nxl 

nxn 

nxl 

mxn 

mxm 

nxm 


TOTAL  STORAGE  *  (M-l)  Mld  +  M2{J  where 


m  -  +  5n 

Mld  "  IT  +  T 


STORAGE  REQUIREMENTS 

2 

n 

Mn 

Mn 

ran 

2 

m 

run 


Mjd  *  j"  +  2 run  +  +  -j— 
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TABLE  9  (continued) 


II.  MULTIPLICATIONS  TO  UPDATE  STORED  MATRICES 


EQUATIONS 

NUMBER  OF  MULTIPLICATIONS* 

1.  Eq. 

(168) 

for  ¥(k) 

n3 

2.  Eq. 

(172) 

/V 

for  a  y j  00 

,  k-M<j  <k 

M  2nm 

— 

3.  Eq. 

(173) 

for  APj(k) 

,  k-M  <  j  <k 

M  2n  m 

4.  Eq. 

(174) 

for  dj (k) , 

k-M<j  dt 

M[nm  +  m] 

5.  Eq. 

(169) 

for  tfi(k) 

2 

n  m 

6.  Eq. 

(170) 

for  A  Vj (k) 

,  k-M<j<k 

2 

M[n  m  +  nm] 

7.  Eq. 

(171) 

for  AK j  (k) 

,  k-M<j  <k 

M[n2  +  n] 

TOTAL  MULTIPLICATIONS  *  (M-2)  N 


2  2 

N.  ,  =*  3n  m  +  n  +  4nm  +n  +  m 

lu 


Id  +  N2d  Where 


N 


2d 


3  2  2 

n  +  7n  m  +  2n  +  8nm  +  2n  +  2m 


*The  sequentially  updated  Kalman  filtering  technique  is  used. 
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Note  the  differences 
Mlc  "  Mld  =  n2 
m2c  "  M2d  =  n2 

Therefore,  the  decomposition  (42)  or  (122)  provides  a  savings 
in  storage  of  (M-l)n2. 

Let  the  number  of  multiplications  be  denoted  by  Nd  for 
the  new  GLR  algorithm  composed  of  a  bank  of  constant  A-state 
filters  using  the  decomposition  of  the  maneuver  signature 
matrix.  The  number  of  multiplications  needed  to  update  the 
stored  matrices  are 

Nd  =  (M-2)  Nld  +  N2<j 


where 


Nid  =  3n2m  +  n2  +  4nm  +  n  +  m 

N2d  =  n^  +  7n2m  +  2n2  +  8nm  +  2n'+  2m 

where  we  have  assumed  that  all  elements  of  the  measurement 
vector  are  uncorrelated  and  the  sequentially  updated  Kalman 
filtering  technique  is  used. 


(183) 

(184) 


(185) 

(186) 
(187) 
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Note  that  Nid  *-s  not  a  function  of  the  third  power  of 
n  and  that 

Nic  -  Nia  *  n3  +  2n2m  (188) 

n2c  ”  N2d  *  2n2m  (189) 

Consequently,  the  decomposition  provides  a  savings  of 

(M-2) [n3  +  2n2m]  +  2n2m  (190) 

for  the  approach  of  using  a  bank  of  constant  A-state  filters. 

Additional  savings  in  multiplications  are  realized  if 
the  process  noise  Q(k)  =0,  for  all  k,  and  if  Y(k)  is  defined  by 

¥(k)  =  Pi(k)  *T(0,k)  (191) 

Note  that  y(k)  appears  only  in  the  defintion  of  AH  (k) ,  Eq. 

(169),  in  the  equations  of  the  constant  A-state  filter.  It 
suffices,  therefore ,  to  compute  AH(k)  without  computing  ¥(k) : 

AH (k )  =  H(k)  Pi(k)  *T(0,k)  (192) 

The  number  of  multiplications  is,  in  general, 

2n2m  (193) 
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If  *(k,0)  is  triangular  that  number  is 

niSnllijT.  (194, 

In  the  above  we  are  assuming  that  no  additional  multiplications 
are  needed  to  compute  $(  0,k)  .  The  matrix  'j'(k)  is  only  needed 
in  (178)  to  obtain  Ax(k)  when  a  jump  has  been  detected.  Since 
jumps  occur  infrequently,  the  matrix  y(k)  needs  computed  infrequently 
for  the  case  of  no  process  noise.  Consequently,  for  the  case 
of  no  process  noise,  the  highest  power  of  n  appearing  in  Table 
9  is  two.  The  total  number  of  multiplications  required  is 
at  most 

5Mn2m  +  2n2m  (195) 

In  this  noiseless  case  the  matrix  YOc)  does  not  need  to  be  stored 

2 

which  results  in  an  additional  storage  savings  of  n  . 
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8.  THE  A-STATE  FORMULATION  FOR  A  PRIORI  JUMPS 

Define  the  A-state  Ax(k,k-1)  as 
Ax(k,k-1)  *  X2(k,k-1)  -  xi(k,k-l),  k^p 

This  random  variable  satisfies 

Ax(k,k-1)  »  ♦(k,k-l)  Ax(k-l),  kj>p 

and,  in  particular 

Ax  (p,p-l)  =  Axp  =  ${q+l,q)  Axq 

Using  (18)  and  (19)  we  see  that  the  state  equation  for  Ax(k,k 
is  given  by 

Ax(k+l,k)  =  A$w  (k+l,k)  Ax(k,k-1),  k>p 

with  initial  condition  (198)  where 

A*w(k+l,k)  =*(k+l,k)  [I  -  Ki(k)  H(k)] 

The  measurement  equation  for  Ax(k,k-1)  satisfies 
Y^(k)  =  H(k)  Ax(k,k-1)  +  YjJk),  k>p 

since 


(196) 

(197) 

(198) 

-1) 

(199) 

(200) 

(201) 


>l(k)  -  Y2(k),  all  k 


(202) 


The  sequence  Yi(k)  is  Gaussian  white  noise,  117],  with  zero 


mean  and  covariance 

E  {Y  i  (k)  YiT(j)}  =  Vi(k)  6kj  all  k 

Note  that 

Ax  (k )  =  [I  -  K1(k)  H  (k )  ]  A  x  (k , k-1) 

A 

Consequently,  if  we  have  the  optimal  estimate  A x (k , k-1 I k) 

A  A 

the  optimal  estimate  Ax{k)  =  Ax(k|k)  is  given  by 

A  A 

Ax  (k)  =  [I  -  Ki(k)  H  (k)  ]  Ax(k,k-l|k) 

Define  the  nxn  matrix  yw(k)  as 

*WU>  =  I 

fwlk)  =  A$w(k,k-1)  y w (k-1)  ,  k>l 

The  matrix  ’J'w(k)  is  positive  definite  for  all  k;  it  has  inverse 
yw“l(k).  We  define  a  new  constant  A -state  variable  A y (k, k-1) 
as 


Ax  (k,k-l)  =  Yw(k)  Ay  (k,Jt-l) 

It  satisfies  the  constant  A -state  equation. 
Ay  (k+l,k)  *  A  y  (k ,  k-1)  ,  k>p 


(203) 

(204) 

(205) 

(206) 
(207) 


(208) 

(209) 
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Its  measurement  equation  is  given  by 
Yl(k)  *  A  Hw  ( k )  Ay  (k , k-1)  +  Yi(k) 

where 

A  Hw ( k )  =AH  (k)  ¥w(k) 

Comparing  (201)  with  (74)  we  find  that  the  maneuver 
signature  matrix  G(k;p)  satisfies 

G(k;p)  A  xp  =  H  (k)  V  w(k)  Vw"l(p)  A  xp 

or,  equivently, 

G(k;p)  Axp  =  Gi(k)  G2<k)  A  xp 

where 

Gi(k)  *  H ( k )  ¥w(k) 

G2(P)  =  ^  w" 1  (P ) 

Note  that 

G2(P)  Ax  -  Ay  (p,p-l) 

Consequently,  Eq.  (213)  may  be  written  as 
G(k;p)  A  x  p  *  A  Hw ( k )  Ay(^*k-1) 


(210) 

(211) 

(212) 

(213) 

(214) 

(215) 

(216) 

(217) 
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Equations  (214)  and  (215)  constitute  a  decomposition  of  the 
maneuver  signature  matrix  into  a  matrix  Gi(k)  which  does  not 
depend  on  the  jump  time  p  and  a  matrix  G2(p)  which  does. 

The  product  of  the  matrix  G2(P)  with  Axp  is  a  constant  as 
shown  by  Eq.  (216) . 


For  the  case  that  the  process  noise  Q(k)  =  0  for  all 
k,  we  can  take  yw(k)  to  be  defined  by 


'Kw(k)  =  Pi(klk-l)  *T(0fk),  k>0 


We  rewrite  (210)  in  the  form 


Azw(k)  =  AHw(k)  Ay(k,k-1)  +  Au(k) 


(218) 


(219) 


1 


where 


Azw(k)  =  Yi (k) 

(220) 

E  (Au(k)  AuT ( j ) }  =  ARwfk) 

(221) 

ARw(k)  =  VX(k) 

(222) 

The  filters  of  the  two  formualtions  (the  a  posteriori  and 
the  a  priori)  differ  only  in  the  defintion  of  the  inputs 
A$  ,  ¥,  Az  and  AR.  Both  formulations  are  equivalent  and 
their  corresponding  A-filters  provide  optimal  estimates  that 


satisfy  the  identify  (205).  The  optimal  estimate  Ax(k,k-1) 
is  obtained  from  AAy(k,k-l)  by  using  (208). 


9.  NORMALIZED  A-MEASUREMENT  NOISE  COVARIANCE 

The  A -measurement  noise  covariance  matrix  AR(k)  (given 
by  (25)  for  the  a  posteriori  jump  formulation  and  by  (222) 
for  the  a  priori  jump  formation)  can  be  normalized  to  R(k) 
by  premultiplying  Eq.  (21)  or  Eq.  (201)  by  the  normalizing 
matrix.  For  this  purpose  we  need  a  matrix  square  root. 

Since  Vi(k)  is  a  symmetric  positive  definite  matrix  it 
may  be  written  in  a  square  root  factored  form  (Kaminski,  et.  al. 
[19])  : 

V(k)  =  V*5  <k) .  V*1  (k)  (223) 

where  V 35  is  a  lower  triangular  matrix  (zeros  above  the  diagonal). 
Square  roots  are  not  necessarily  unique  but  a  unique  square 
root  may  be  defined  using  the  Cholesky  decomposition.  We 
assume  the  Cholesky  decomposition  is  used  and  denote  the  square 
root  of  a  matrix  with  the  superscript  ^  and  its  transpose 
by  %t.  The  square  root  of  the  inverse  of  V(k)  is  denoted  by 
V_Js(k)  . 

Premultiplying  (21)  by  the  factor  R ^(kjV*57  (k)  R-1(k) 
and  redefining  Az(k),  AH(k)  and  Au  (k)  as 

Az  (k)  =  R*(k)V*T(k)  R_1(k)  [z(k)  -  H(k)  xi(k)]  (224) 


61 
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AH (k )  =  R35(k)vJ4T(k)  R_1(k)  [H (k )  J 

Au  (k )  =  R3s(k)  v^U)  R-1  (k)  [z(k)  "  H(k)  x2<k)] 


(225) 

(226) 


we  obtain  the  form  (21)  where  the  new  Au(k)  is  a  Gaussian 
white  noise  sequence  with  zero  mean  and  covariance 


E  {Au (k)  AuT(j) }  =  AR(k)  6kj 


(227) 


where 


AR(k)  =  R(k) 


(228) 


which  is  the  covariance  of  the  measurement  noise  v (k)  of  (5) . 


The  normalizing  factor  for  the  a  priori  jump  formulation 
(201)  is 


R^k)  Vi_J5(k) 


(229) 


The  normalization  is  particularly  useful  for  the  case 
that  R(k)  is  a  constant.  If  there  is  no  jump,  Az(k)  is  zero 
mean  with  covariance  R(k) .  If  there  is  a  jump,  Az(k)  has 
mean  AH(k)  Ax(k)  and  covariance  R(k).  Looking  for  a  jump 
may  be  avoided  if  the  residuals  Az(k)  appear  to  be  zero  mean 
with  covariance.  R()c)  . 
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10.  EXTENSION  OF  APPROACH  TO  CONTINUOUS  LINEAR  STOCHASTIC 
SYSTEMS 


Consider  the  continuous  linear  discrete  stochastic  system 

A 

described  by  the  vector  (Ito)  stochastic  differential  equation 

dx(t)  =  F(t)  x ( t )  dt  +  T (t)  dw(t),  t>  0  (230) 

Ak 

and  the  vector  (Ito)  observation  equation 

dZ  (t)  =  H  (t)  x(t)  dt  +  du(t)  ,  t  >  0  (231) 

where  x(t)  is  the  n-vector  state,  F(t)  and  r  (t)  are,  respectively, 
nxn  and  nxr  nonrandom,  continuous  matrix  time-function,  and 
(w(t),  t>0 }  is  an  r-vector  Brownian  motion  (Wiener)  process 
with 

E{ dw ( t )  dw(t)T}  =  Q ( t )  dt  (232) 

The  observed  process  { Z ( t ) ,  to  0}  is  an  m-vector  process,  H(t) 
is  an  mxn,  nonrandom,  continuous  matrix  time-function,  and 
{u (t) ,  t>0}  is  an  m-vector  Brownian  motion  (Wiener)  process 
with 

E{dy(t)  du(t)T}  =  R(t)  dt  (233)' 

where  R(t)>0  .  We  assume  that  the  system  (230)  and  (231)  is 
observable. 
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The  well-known  continuous  Kalman-Bucy  filter  [14,15] 


is  given  by,  [ 16] , 


x  ( t )  =  F  (t)  x  ( t )  +  K  (t)  [ 2  ( t )  -  H  ( t)  x  ( t )  ]  ,  t>  0 


P 


(t) 


=  f (t)  p(t)  +  p(t)  F?(t)  +  r(t) 

-  K  ( t) 


Q (t)  rT(t) 

H  ( t )  P(t)  ,  t>  0 


where 


K ( t )  =  P(t)  HT(t)  R-l{t) 

and  where  the  formalized  dZ(t)/dt  is  written  as  z(t).  We 

A.  A 

have  made  the  identifications  x(t)  =  x(t|t)  and  P(t)  =  P(t|t). 

A  jump  in  state  with  magnitude  £Xq  occurs  at  time  g: 
x(t+)  =  x (t)  +  AXq  for  t  =  q 

Consider  the  three  Conditions  Hi,  Hi  and  H2-  Let  both  x(t) 
and  X2(t)  represent  the  state  for  the  system  described  by 
(230),  (231)  and  (237).  Let  xi(t)  represent  the  state  for 
the  case  that  the  jump  magnitude  AXq  is  zero  (i.e.,  there 

*  ~  /v 

is  no  jump).  Let  xi(t),  xi(t)  and  X2(t)  denote  the  Kalman-Bucy 
filter  estimates  of  the  state  (i.e.,  Eqs.  (234)  -  (236))  for 
the  Conditions  Hi,  Hi  and  H2,  respectively.  The  relationship 
between  these  estimates  are  similar  to  those  given  in  Figure  1 
for  the  discrete  system.  The  gains  for  the  three  Conditions 
are  equal  as  well  as  the  covariances. 


(234) 


(235) 


(236) 


(237) 
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The  estimate  xi(t)  can  be  interpreted  as  follows.  It 

a 

is  identified  with  X2(t)  up  to  time  q 
xi (t)  =  X2 ( t)  t<q 

At  time  q+  the  amount  AXq  is  subtracted  from  the  filtered 
estimate  X2(t)  and  the  result  is  defined  as 

A 

xi(t)  =  X2(t)  -  ax q,  t  *  q+ 

For  t>q  the  quantity  satisfies  Eq.  (234). 


We  are  interested  in  the  difference  between  X2(t)  and 
xi (t)  for  t>q.  We  define  this  difference  as 

Ax(t)  «  X2(t)  -  xi(t),  t>q 

It  has  the  initial  condition 
AX(t)  =  Axq,  t  =  q+ 


Since  both  X2<t)  and  xi(t)  satisfy  (234),  the  difference 
Ax(x.)  satisfies 


dAx  (t) 
*t 


tF(t)  -  Ki(t)  H  (t)  ]  Ax  (t)  ,  t>q+ 


(238) 


where  we  have  used  the  identity 
22(fc)  *  Zl(t) 


(239) 
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Define  the  nxn  nonrandom,  continuous  matrix  time-function 
Y(t)  as  the  solution  to  the  differential  equation 

¥<t)  =  [F  ( t)  -  Ki(t)  H  (t)  ]  y(t),  all  t 

Y(  0)  =  I 

The  matrix  time-function  y(t)  is  positive  definite  and  has 
inverse  y"l(t).  Define  the  new  A~state  variable  Ay(t)  as 

AX(t)  *  Y(t)  Ay(t),  t>q+ 

The  substitution  of  (242)  into  (238)  gives 

^  ■  0 . tv 

We  note  that  &y  (t)  is  a  constant  and  satisfies 
AY  (t)  =  '?-1(q)  AXq,  t>q+ 

The  measurement  equations  for  a*  and  AY  are  given  by 
Az(t)  =  H ( t )  Ax(t)  +  Au(t) ,  t>q+ 

Az(t)  *  AH (t)  A  Y(t)  +  Au(t) ,  t >q+ 


where 


AH  (t)  =  H  (t)  y(t)  ,  all  t 


Az(t)  «  z(t)  -  H(t)  xi  (t)  ,  all  t 


* 

Au(t)  -  Yi(t)  -  Z1(t)  -  H ( t)  xi(t),  all  t 

«Y2(t)  -  z2(t)  -  H (t)  x2(t),  all  t 

an 


(240) 

(241) 

(242) 

(243) 

(244) 

(245) 

(246) 

(247) 

(248) 


(249) 


Eqs.  (247)  and  (244)  constitute  the  decomposition  of  the  maneuver 
signature  matrix.  The  residual  Yi(t)  is  the  measurement  noise 
for  the  A~process.  It  is  zero  mean  and  it  has  the  same  covariance 
as  u(t) ,  117] : 

E{Au(t)  AuT(s)  }  -  a  R(t)  6  (t-s)  (250) 

where  fi(t-s)  is  the  Dirac  delta  function  and 

AP(t)  =  R(t)  (251) 

If  the  jump  time  q  were  known  and  if  the  initial  state 

A 

AY(q)  were  normally  distributed  with  mean  AY(q)  and  covariance 
AP(q)  then  the  Kalman-Bucy  filter  applied  to  (243)  and  (246) 
would  provide  the  solution: 

m 

AY(t)  =  AK(t)  [Az(t)  -  AH  (t)  AY(t)]  (252) 

AP(t)  =  -A  K  ( t )  AH  (t)  AP(t)  (253) 


where 

A  K(t)  =  A  P  ( t)  AHT(t)  R'ijt)  (254) 

A 

Using  (242)  the  optimal  estimate  x(t)  of  x(t)  is 

x  ( t )  *  xj_  ( t )  +  y  (t)  Ay(t)  (255) 

with  covariance 

P(t)  *  Pi(t)  +.  f(t)  A  P (t)  Y  T ( t)  (256) 
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The  jump  time  q  is  unknown  and,  therefore,  must  be  estimated 
along  with  Ay.  For  each  t>q,  we  desire  to  obtain  the  estimates 

A  A  A 

q(t)  and  Ay  (q(t),t)  that  render  a  minimum  to  the  function 

q 

J(q»  Ay;t)  *  /  Az(t)t  AR(x)-1  Az(t)  dx 

0 

t 

+  /  (Az  (t)  -  AH (t )  Ay]T  AR(t)-1  (Az(x)  -  AH (x  )Ay]  dr 

q 

(257) 


or,  equivalently,  a  maximum  to  the  log  likelihood  ratio 

t 

Mq,Ay;t)  *  /  Az(t)t  AR(t)"1  Az(t)  dx  -J(q,Ay;t) 

0 

(258) 


The  above  integrals  are  Ito  integrals.  For  a  fixed  q,  the 
optimizing  Ay(q,t)  satisfies 


C (q; t)  A  y (q, t)  =  D (q; t) 


where 


t 

C(q;t)  *  /  AH(x)T  AR(x)_1  AH(x)  dx 

q 


D(q;t)  =  fl  AH(x)t  AR(x)-1  Az(x)  dx 

q 


(259) 

(260  ) 

(261) 
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We  make  the  definitions 


c  ( t )  *  AH(t)T  dR(t)-l  A  H  ( t)  (262) 

d (t)  *  4H(t)T  A  R(t)“l  A  z(t)  (263) 


we  note  that  d(t)  is  an  nxl-vector  Brownian  motion  process. 


In  view  of  (259),  Eq.  (258)  becomes 


«.  (q»  A  y  (q»t)  ;t)  =  DT(q;t)  C_1(q?t)  D  (q;t) 


(264) 


=  I  /  dT(T)  dT)  C-l(q;t)  [ 

q 


f  d(x)  dx  ] 

q 

(265) 


The  maximizing  argument  of  (265)  is  denoted  by  q(t) . 

Returning  to  Eqs.  (255)  and  (26)  the  Kalman  gain  K(t) 
satisfies 

K(t)  =  Kx(t)  +  y  (t)  A  K ( t) 

The  continuous  case  addressed  in  this  section  is  being 
treated  in  more  detail  in  another  report. 
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11.  THE  SIMILARITY  OP  JUMP  ESTIMATION  AND  BIAS  ESTIMATION 


The  bias  estimation  problem  can  be  expressed  as  follows, 

[20], 


State'  Dynamics 

x  *  Ax  +  Bb  +  w  (266) 

Bias  Dynamics 

b  *  0  (267) 

Observation  Equation 

z  «  Hx  +  Cb  +  u  (268) 

where  the  state  x  is  an  n-vector,  the  bias  b  is  a  p-vector, 
the  observation  vector  y  is  an  m-vector,  w  is  the  process 
noise  vector  with 

E [w ( t )  wT(s)]  =  Q ( t)  6 (t-s) 

and  u  is  the  observation  noise  vector  with 
E[u(t)  uT(s)]  *  R ( t )  6  (t-s) 

The  vectors  w  and  u  are  assumed  to  be  independent.  The  matrices 
A  and  B  are  time  varying. 
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Friedland  [20]  was  the  first  to  show  that  the  optimal 

A 

estimate  x(t)  of  the  state  x(t)  could  he  expressed  as 

x(t)  *  X1(t)  +  S(t)  b  (t)  (269) 

where  X].(t)  is  the  bias-free  estimate  '  i  •  s . ,  the  filter  assumes 
b  is  zero  even  though  it  is  nonzero),  h(t)  is  the  bias  estimate 
which  is  computed  using  the  bias-free  residuals  z(t)  -  H(t)  xi(t) 
and  S (t)  satisfies  the  differential  equation 

S(t)  =  [A (t)  -  Ki(t)  H (t) ]  S (t)  +  B (t)  -  Ki(t)  C ( t)  (270) 

with  initial  condition 

S(0)  =  0  (270) 

The  matrix  Ki(t)  is  the  Kalman  gain  for  the  bias-free  estimate. 

A 

The  covariance  P(t)  of  x(t)  is  shown  in  [20]  to  satisfy 

P(t)  =  Px(t)  +  S  (t)  Pb(t)  ST(t)  (271) 

where  Pi(t)  is  the  covariance  of  xj.(t)  and  Pb(t)  is  the  covariance 
of  b. 


Eqs.  (255)  and  (256)  have  the  same  form  as  (269)  and 
(271).  The  bias  estimation  problem  is  equivalent  to  the  jump 
estimation  problem  of  estimating  the  jump  b(0)  which  occurs 
at  time  0;  that  is,  the  value  of  the  jump  state  b  is  zero 
before  the  jump  ajid  b(0)  afterwards.  The  jump  time  0  is  known. 
The  bias  estimation  problem  for  discrete  systems  is  also  treated 

in  [20].  An  extension  to  the  problem  of  indirect  observations 
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is  given  in  [24] . 


Friedland's  bias  filtering  technique  is  extended  in  [21] 
to  the  case  of  estimating  a  time  varying  bias 

b  =  F (t) b  (272) 

The  solutions  have  the  forms  (269)  and  (271) . 

Mendel  and  Washburn  [22]  show  that  the  estimation  of 
the  bias  vector  b  can  be  interpreted  as  being  equivalent  to 
the  estimation  of  a  constant  that  is  observed  through  white 
noise.  That  result  compares  with  Eqs.  (43)  and  (246).  That 
interpretation  is  reviewed  in  [23] . 

The  structure  (269)  is  shown  in  [25]  to  hold  for  the 
optimal  state  estimate  under  the  uncertainty  of  different 
failure  modes. 

The  problems  of  jump  estimation  and  bias  estimation  differ 
in  the  following  way.  The  bias  estimation  is  that  of  estimating 
a  parameter  which  undergoes  a  single  jump  from  a  zero  value. 

The  jump  estimation  problem  is  that  of  estimating  a  parameter 
that  undergoes  multiple  jumps  from  nonzero  values.  Once  a 
jump  has  been  detected  and  estimated  it  is  necessary  for  the 
jump  estimator  to  pass  this  information  on  to  the  original 
state  estimator  and  to  reinitialize  for  the  next  jump.  The 
estimator  described  by  (269)  is  for  single  jump  systems. 

The  estimator  described  by  (255)  is  for  multiple  jump  systems. 


A  nonlinear  algorithm  is  given  in  [26]  -  [28]  for  detecting 
and  estimating  sudden  changes  of  biases  in  linear  stochastic 
systems.  The  method  is  based  on  maximum-likelihood  estimation 
[29]  . 

An  extension  of  Friedland's  bias  filtering  technique 
to  a  class  of  nonlinear  systems  is  given  in  [30). 


12.  CONCLUSION 

\ 

\ 

We  have  presented  two  recursive  GLR  algorithms  for  detecting 
and  estimating  maneuver  states  and  parameters  in  the  engagement 
between  an  anti-ship  cruise  missile  and  a  ship  defense  interceptor. 
A  decomposition  of  the  maneuver  signature  matrix  is  used  to 
derive  the  algorithms.  The  computational  and  storage  requirements 
are  substantially  less  than  those  of  other  GLR  algorithms.  The 
decomposition  divides  the  maneuver  signature  matrix  into  the 
product  of  two  matrices.  One  matrix  depends  only  on  the  current 
observation  time  while  the  other  depends  only  on  the  jump  time. 

The  product  of  the  latter  matrix  and  the  jump  magnitude  vector 
provides  a  jump  error  state  vector  which  is  constant.  This 
constancy  facilitates  using  the  GLR  approach.  The  other  matrix 
of  the  decomposition  represents  the  new  maneuver  signature  matrix 
for  the  new  constant  jump  error  state  vector.  The  nondependency 
of  the  maneuver  signature  matrix  on  the  jump  time  avoids  storing 
large  matrices  and  computing  large  matrix  products  for  each 
past  observation  time. ^ 

Previously  the  most  efficient  GLR  algorithm  for  the  detection 
and  the  estimation  of  jumps  required  (at  each  current  observation 
time)  multiplications  on  the  order  of 

3  2 

Mn  +  5Mn  m 
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where  n  is  the  dimension  of  the  state  vector,  m  is  the  dimension 
of  the  measurement  vector  and  M  is  the  number  of  candidate  jump 
times  in  the  past.  As  a  consequence  of  the  decomposition  the 
GLR  algorithm  II  presented  herein  requires  multiplications  on 
the  order  of 

n^  +  3Mn^m 

for  general  discrete  linear  stochastic  systems  and 
5Mn^n 

for  such  systems  without  process  noise.  This  is  a  substantial 

savings  in  computation.  The  corresponding  savings  in  storage 

2 

is  on  the  order  of  Mn  . 

Algorithm  II  is  a  bank  of  Kalman-Bucy  constant  A -state 
filters  that  use  the  decomposition  of  the  maneuver  signature 
matrix.  There  is  a  filter  for  each  candidate  jump  time  in  the 
past.  As  a  result,  the  required  computations  for  jump  detection 
and  estimation  are  performed  at  each  observation  time.  In  contrast, 
algorithm  I,  using  minimal  computations,  computes  and  stores 
at  the  current  observation  time  an  information  matrix  for  later 
processing.  The  appropriate  information  matrix  for  any  candidate 
jump  time  is  obtained  simply  by  taking  the  difference  between 
the  information  matrix  of  the  current  time  and  that  of  the  candidate 


jump  time.  The  required  computations  for  jump  detection  and 
estimation  are  performed  in  the  vacinity  of  a  jump  or  false 
alarm.  While  algorithm  I  requires  more  computations  at  an  observat 
time  than  algorithm  II  to  detect  and  estimate  jumps  it  is  not 
necessary  that  the  calculations  be  performed  at  each  observation 
time.  Consequently,  algorithm  II  is  more  applicable  for  systems 
with  frequent  jumps  and  algorithm  I  more  applicable  for  systems 
with  infrequent  jumps  or  maneuvers. 

The  GLR  algorithms  I  and  II  are  a  computation  improvement 
over  existing  techniques,  algorithms  and  methods  [1,2],  [4-6], 

[8],  [9],  and  [20-38]  for  adaptively  estimating  the  state  of 
a  linear  stochastic  system  undergoing  abrupt  changes  (e.g., 
maneuvers)  in  state. 
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APPENDIX  A 


RELATION  BETWEEN  KALMAN  GAINS 


The  argument  of  the  following  variables  is  k:  1 
Vj,  AV,  K,  Kj  and  AK.  The  argument  of  A  P  is  k-1  and 
of  *  is  (k,k-l) .  The  argument  of  P  and  is  (k|k-l) 

we  have  from  Eq.  (B-7) 

V  *  Vx  +  H  *  AP  $T  HT 

Postmultiplying  (A-l)  by  V-1  gives 

T  T  -1  ~  -1 

H*AP$aH  V«I-  VXV 

Premultiplying  (A-2)  by  Kj  gives 
KjH^AP^Vv"1  *  K1II-V1V"1] 

The  matrix  AK  satisfies 

akav  «  [i-k1h]*ap«thtv1“1r 

-1  ~  -1 

Postmultiplying  (A-4)  by  AV  RVj  gives 

~  -1  T  T  -1  T  T  -1 

AKRV1  1  -*AP*  H  V  i-K1H*AP*  H  V  1 


A-l 


,  R»  V, 
the  argument 


(A-l) 


(A-2) 


(A- 3) 


(A-4) 


(A-5) 


Substitution  of  (A-3)  into  (A-5)  gives,  in  view  of  (B-9) , 


AKII-HKJ  »*AP*THTV“1  +  K1V1V~1-K1  (A-6) 

Since 

K  =  PHTV_1  (A-7) 

and  since,  from  (B-6) , 

T  -1  T  -1  T  T  -1 

PHaV  a  =  P1HiV  A  +4APJ  rV  1  (A- 8) 

it  follows  that 

K  =  K^V-1  +$AP*THTV~1  (A- 9) 

Substituting  (A-9)  into  (A-6)  gives 

K  -  Kx  +  AKfl-HKj]  (A-10) 

A-2 


APPENDIX  B 


RELATION  BETWEEN  PREDICTED  MEASUREMENT  RESIDUAL  COVARIANCES 

The  matrices  V(k),  (k)  ,  P(k|k-1)  and  P1(k|k-1)  are  given  by 

V(k)  =  H(k)  P  (k|  k-1)  HT  (k)  +R(k)  (B-l) 

V^k)  =  H  (k)  Px  (k|  k-1)  HT(k)  +  R(k)  (B-2) 

P(k|k-1)=  4>(k,k-l)  P  (k-1)  4  T(k,k-1)  +T  (k-1)  Q  (k-1)  T  T  (k-1) 

(B-3) 

Px (k | k-1)  =$(k,k-l)  P1(k-1) *  T(kfk-1)  +  T(k-l)  Q(k-l)  rT(k-l) 

(B-4 ) 


Substituting 

P (k-1)  =  P1(k-1)  +  AP (k-1)  (B— 5 ) 

into  (B-3)  gives,  in  view  of  (B-4) , 

P (k | k— 1 )  *  Px  (k j k-1)  +  *(fc,k-l)  AP(k-l)  *T(k,k-l)  (B-6) 

Substituting  (B-6)  into  (B-l)  gives,  in  view  of  (B-2), 

V  (k)  *  Vx(k)  +  H(k)$  (k,k-l)  AP  (k-1)  *T(k,k-l)  HT(k)  (B-7) 

The  matrix  AP(k|k-l)  is  given  by 

A P (k| k-1)  -  [i-K^kjHWH^k-DAPfk-l)*  T(k,k-1) 

[MjtWHlk)]1  (B-8) 


B-l 


\ 

\ 

\ 


Recall  the  identity 

R(k)V1“1(k)  =  tl-H(k)  Kx (k) ]  (B-9) 

Premultiplying  (B-8)  by  H(k)  and  postroultiplying  that  product  by 
HT(k)  gives,  in  view  of  (B-9), 

H(k)  A  P (k| k-1)  HT (k)  =  R(k)V1"1(k)  F (k,k-l) V1 (k)  R(k)  (B-10) 

where 

F  (k,k-l)  =  H  (k)4>  (K,k-1)  AP(k-l)  4>T(k,k-l)  HT(k)  (B-ll) 

The  matrix  A V(k)  is  given  by 

A  V (k)  =  H(k)  AP  (k|  k-1)  HT(k)  +  R(k)V1"1(k)  R(k)  (B-12) ' 

Using  (B-12)  to  ~'lve  for  F(k,k-1)  in  (B-10)  gives 

F  (k,k-l)  *V1(k)R"1(k)  AV(k)  R_1  (k)  V1  (k)  -  V1  (k)  (B-13) 

Equating  (B-ll)  and  (B-13)  and  substituting  the  result  into 
(B-7)  gives 

V(k)  *V1(k)  R-1  (k)  Av(k)  R_1  (k)  Vx(k) 


(B-14 ) 


APPENDIX  C 


RELATION  BETWEEN  STATE  ESTIMATES 

We  desire  to  show  that  the  state  estimates  satisfy 

A  ^  A 

x(k)  =  x1 (k)  +  Ax(k)  (C-l ) 

We  assume  that 

x(k-l)  =  x1(k-l)  +  Ax(k-l)  (C-2) 

The  individual  estimates  satisfy  the  following  equations 

x (k)  =  [I-K(k)H(k)l  $  (k,k-l)  x(k-l)  =  K(k)z(k) 
xx(k)  =  [I-K^ (k) H (k) ]  $ (k,k-l)  x1(k-l)  +K1(k)z(k) 

Ax (k)  *  [I-AK(k)H(k)]At  (k,k-l)  A*(k-1)  +  a  K (k) 

[z(k)  -  H (k)  xx (k) ] 

Using  (C-2)  -  (C-5)  it  follows  that  (C-l)  is  satisfied  provided 

K (k)  »  K1  (k)  +  A K (k)  [I-H(k)  Kj  (k)  ]  (C-6) 

(I-K(k)H(k) ]  $  (k,k-l)  =  [I-  a K (k)H (k) ]  a*  (k,k-l)  (C-7) 

Substituting  the  definition  of  A$  (k,k-l)  into  (C-7)  results  in  an 
equation  which  holds  provided  (C-6)  holds.  The  validity  of  (C-6) 
follows  from  Appendix  A. 


(C-3) 

(C-4) 

(C-5) 


I 

1 

APPENDIX  D 

RELATION  BETWEEN  COVARIANCES 


We  wish  to  show  that 
P  (k)  =  P1  (k)  +  AP(k) 

given  that 

P(k-l)  =  Px(k-1)  +  AP(k-l) 

The  argument  of  t:he  matrices  H,  K,  ,  AK  and  AR  is  k.  The 
argument  of  $  and  A$  is  (k,k-l) . 

The  covariance  P (k)  satisfies 
P (k)  =  [I— KH ]  P (k| k-1) 

From  Eq.  (C-7)  we  have 

[I-KH]  =  [I-AKH]  [I-KjHJ 

From  Eq.  (B-6)  we  have 

P (k| k-1)  =  P1(kjk-1)  +  $AP  (k-l)$T 


(D-l) 


(D-2) 


(D-3) 


(D-4) 


(D-5) 


1 

I 

! 


D-l 


~T 


*BRos;o55^^ISr^2'0’ 

JUL  81  ML  5>ALF0RD  ^ 

UNCLASSIFIED  PSI-TR-S1-1  _ 


(D-6) 


The  covariance  P1 (k)  satisfies 
P1(k)  -  [I-^HJ  P1(k|k-1) 

Substitution  of  (D-4)  -  (D-6)  into  (D-3)  and  making  use  of  the 
definition 

M  =  [I-K1H]  $  (D-7) 

results  in 

P (k)  -  [I  -  AKH]  [Pj (k)  +  A*AP(k-l)*  T)  (D-8) 

or ,  equivalently , 

P  (k)  *  P1(k)  -  AKHP1  (k)  +  H-AKH]  A4>AP  (k-1)  ♦  T  (D-9) 

The  covariance  AP(k)  satisfies 

AP (k)  =  [I- AKH]  A*AP(k-l)A*  T  (D-10) 

Making  use  of  (D-7) ,  Eq.  (D-10)  can  be  written  as 

[I-AKH]  A*AP(k-l)*T  ■  AP(k)  +  AP  (k)HTR_1V1K1T  (D-ll) 


1 

D-2  I 


since 


V1"1  R  *  [I-HKj^  (D-12) 

It  follows  from  (D-12)  and 

H  P2(k)  -  R  KxT  (D-13) 

That 

AKHP1(k)  *AK  AR  R (D-14) 
where 

Ar  «  R  V1~1  R  (D-15) 

Substitution  of  (D-ll)  and  (D-14)  into  (D-9)  gives 

P(k)  =P1(k)  +  A  P  (k)  +  [A  P  (k)  HT  -  AKARJR-^^1  (D-16) 


Eq.  (D-l)  now  follows  since 

AP(k)HT  *  AK  AR  (D-17) 


D-3 


APPENDIX  E 


CHANG  AND  DUNN'S  SYSTEM  FORMULATION: 
DISCUSSION  OF  IMPULSE  INPUT  VERSION 


The  following  system  step  input  case  is  considered  in  [8] 


and  [9]: 

X (k+1)  *  A(k+l,k)  X (k)  +  B (k+1 ,k) u (k,q)  +  r(k)w(k)  (E-l) 
u (k+1 ,q)  =  F (k+1 ,k)  u(k,q)  (E-2) 
u(k,q)  *  0  for  k<q  (E-3) 
u(q,q)  #  0  (E-4) 


where  F  is  the  state  transition  matrix  for  the  control  variable 
u  and  the  other  variables  are  as  defined  in  Eqs.  (1)  and  (3)  of 
Section  1.  The  unknowns  to  be  estimated  and  detected  are  u(k,q) 
and  q. 

We  note  two  differences  between  the  dynamics  defined  by 
(1)  and  (2)  and  that  defined  by  (E-l)  -  (E-4) .  First,  the  matrix 
F  is  taken  as  the  identity  matrix  in  Eq.  (2) .  Secondly,  the  value 
of  u  is  zero  before  the  jump  in  the  above  formulation  but  may  be 
nonzero  in  Eq.  (1) . 

Eqs.  (E-l)  -  (E-4)  and  Eq.  (3)  is  the  "step  input  to  the 
ftystem"  case.  It  is  described  in  [3]  as  the  "dynamic  step"  case 
in  which  the  input  matrix  B  and  the  matrix  F  are  identity  matrices. 


The  step  input  case  can  be  transformed  into  the  "impulse 
input  to  the  system"  case  by  augmenting  the  state  vector  X  with 
the  control  vector  u.  Let  x  denote  the  augmented  state  vector. 
The  augmented  system's  dynamics  are  given  by 


x(k+l)  *  *(k+l,k)  [x(k)  +  Axq  6qk)  +  T  (k)  w(k) 


(E-5) 


where 


*  (k+l,k)  »  [ 


A(k+l,k) 


B (k+1 ,k) 
F(k+l,k) 


(E-6) 


ixq  *  [j]  u  (q  ,q) 


(E-7) 


and  where  r  (k)  is  redefined  as  the  augmented  matrix  as  in 
Eq.  (10). 


In  this  report  we  are  interested  in  the  general  formulation 
of  (E-5)  in  which  $,  Axq  and  r  are  arbitrary,  not  necessarily 
satisfying  Eqs.  (E-6) ,  (E-7)  and  (10).  We  desire  to  apply 
Chang  and  Dunn's  GLR  algorithm  in  our  context.  Their  algorithm 
consists  of  a  bank  of  Kalman-Bucy  filters  to  estimate  the  jump 
Axq.  The  filters  are  defined  by  Eqs.  (123)  -  (127) ,  (132)  -  (138) 

A 

and  (142)  -  (143)  in  Section  6.  The  estimate  Avq(k)  represents 
the  optimal  estimate  of  Axq  given  the  measurements  up  to  time  k. 
Consequently,  Chang  and  Dunn's  GLR  algorithm  as  made  use  of  in 
this  report  is  the  impulse  input  version  of  their  technique  as 
described  in  [8]  and  (9]  for  the  "step  input"  case. 


E-2 


For  a  particular  problem  in  the  form  of  (E-l)  -  (E-4)  it 
is  more  efficient  to  use  the  "step  input"  case  filters  as 
described  in  [8]  and  [9]  than  it  is  to  use  the  filters  of  the 
augmented  "impulse  input"  case.  In  this  report,  however,  we 
are  interested  in  the  "impulse  input"  case  when  it  is  not 
necessarily  reducible  to  the  "step  input"  case. 


